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INTRODUCTION 

This  is  the  1979  annual  report  of  the  NRL  Inertial  Confinement  Fusion  Theory  Program.  It  cov¬ 
ers  research  performed  from  October  1978  through  December  1979.  Research  in  each  of  the  four 
current  program  areas  is  reported: 

Section  A.  Laser  Light  Absorption. 

Section  B.  Fluid  Dynamics  of  Ablative  Acceleration 

Section  C.  Development  of  Computational  Techniques. 

Section  D.  Rayleigh-Taylor  Stabilization  Techniques. 

Our  work  in  Section  A  includes  an  analysis  of  anomalous  absorption,  backscatter,  and  flux  limitation  in 
laser-produced  plasmas.  The  physical  processes  included  in  the  model  are  inverse  bremsstrahlung, 
resonant  absorption,  absorption  by  ion-acoustic  fluctuations,  and  Brillouin  backscatter.  In  Section  B  we 
report  our  work  on  the  energy  conversion  of  light  ion  beams  through  the  ablative  acceleration  of  thin 
foils,  our  work  on  spherical  shock  implosion  dynamics  and  stability,  and  our  ID  and  2D  work  on  the 
dynamics  of  laser-driven  thin  plastic  foils.  In  Section  C  we  report  our  current  developments  in  triangu¬ 
lar  grid  techniques.  These  results  include  improved  calculations  of  the  Rayleigh-Taylor  instability 
through  higher  resolution  in  the  SPLISH  numerical  model.  In  Section  D,  we  report  our  improved 
understanding  of  the  "quasi-static"  equilibrium  for  laser-driven  ablation  layers,  and  discuss  the  evolu¬ 
tionary  trends  for  the  various  possible  "quasi-static"  equilibria.  We  also  report  on  our  development  of  a 
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vorticity  generation  model  to  study  the  Rayleigh-Taylor  instability  and  include  our  results  for  test  prob¬ 
lem  and  actual  laser-driven  ablation  layer  profiles.  We  also  include  our  work  on  the  Rayleigh-Taylor 
eigenvalue  problem. 

The  report  concludes  with  a  small  section  on  our  plans  for  the  extension  and  continuation  of  this 
research.  Footnotes,  figures,  and  references  for  each  section  below  are  compiled  separately. 

A.  LASER  LIGHT  ABSORPTION 

The  absorption  and  backscatter  of  light  in  the  underdense  blowoff  of  a  laser  produced  plasma,  as 
well  as  limitiation  on  electron  thermal  energy  flux  all  play  important  roles  in  laser  fusion.  This  year’s 
work  has  been  a  continuation  of  our  earlier  work  in  this  area1,2  and  models  the  dynamics  of  the  interac¬ 
tion  rather  than  the  steady  state.  Also  included  in  this  effort  has  been  a  study  of  the  effects  of  Brillouin 
backscatter.  Although  our  work  applies  in  principle  to  any  laser  produced  plasma,  it  is  particularly 
motivated  by  experiments  at  the  Naval  Research  Laboratory  on  planar  targets.1-6  Experiments  here 
have  demonstrated  that  backscatter  through  the  lens  (apparently  Brillouin  backscatter)  can  be  a  very 
important  process  for  structured  laser  pulses  at  sufficiently  high  intensity.3-4  Also,  more  recent  experi¬ 
ments  with  longer  pulses5  (pulse  duration  of  about  3  nsec)  have  shown  very  high  absorption  even  at 
intensities  where  inverse  bremsstrahlung  is  not  effective.  For  instance  at  an  irradiance  of  7  x  1014 
W/cm2,  the  fractional  absorption  is  about  60%.  Experiments  at  NRL6  and  elsewhere7  also  have  shown 
that  absorption  is  enhanced  when  the  target  is  in  the  focal  plane  of  the  lens.  It  has  always  been  our 
opinion  that  resonant  absorption  alone  cannot  account  for  all  of  the  absorption.2  Reference  2  reviews 
absorption  measurements  done  at  many  laboratories  and  comes  to  this  conclusion.  The  additional 
mechanism  which  we  propose  is  ion  acoustic  turbulence  driven  by  the  return  current.  That  mechanism 
takes  place  in  either  an  unmagnetized1  or  magnetized2  plasma.  The  fluid  simulations  reported  here 
attempt  to  model  this  process.  A  very  much  oversimplified  picture  is  that  ion  acoustic  turbulence  in 
the  underdense  plasma  begins  to  take  over  where  classical  inverse  bremsstrahlung  becomes  ineffective. 
For  instance,  the  solid  line  in  Figure  A.l  is  a  schematic  of  the  total  electron-ion  momentum  exchange 
collision  frequency  as  a  function  of  laser  irradiance  for  a  Nd  laser  produced  plasma.  The  dotted  line  is  a 
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plot  of  the  classical  inverse  bremsstrahlung,  while  the  dashed  line  is  a  plot  of  the  effect  of  the  ion 
acoustic  instability.  Roughly  speaking,  ion  acoustic  turbulence  at  high  temperature  gives  perhaps  half 
as  much  absorption  as  inverse  bremsstrahlung  at  low  temperature. 

One  of  the  principal  conclusion  of  this  work  is  that  the  absorption  physics  falls  into  two  different 
regimes.  In  short  single  pulse  experiments,  resonant  absorption  is  the  dominant  process,  but  it  gets  a 
strong  boost  from  the  ion  acoustic  turbulence.  Typically  the  total  absorption  is  about  40%  with 
resonant  absorption  accounting  for  25%  and  ion  acoustic  turbulence  for  15%.  Most  of  the  light  is  spec¬ 
ularly  reflected,  although  stimulated  Brillouin  backscatter  plays  a  role  also.  On  the  other  hand,  for  long 
pulses  or  structured  pulses,  the  physics  seems  to  be  dominated  by  a  relatively  complicated  interplay 
between  backscatter  and  ion  acoustic  turbulence.  Both  of  these  processes  depend  upon  the  quasi-linear 
and  nonlinear  behavior  of  a  plasma  instability,  so  the  physics  is  necessarily  more  speculative  than,  for 
instance,  absorption  by  inverse  bremsstrahlung.  This  modeling  is  more  complicated  in  a  fluid  simula¬ 
tion  but  can  nevertheless  bring  a  more  complete  picture  of  absorption  in  laser-plasma  interactions.  Our 
treatment  of  these  processes  is  internally  self-consistent  and  is  based  upon  well  documented  theories, 
simulations  and  experiments. 

We  will  now  briefly  contrast  our  approach  to  the  physics  of  the  underdense  plasma  with  particle 
simulations  and  other  fluid  simulations.  Particle  simulations  have  been  particularly  useful  in  studying 
resonant  absorption*- 10,  strongly  driven  backscatter  instabilities  in  a  laser  produced  plasma11  and 
current  driven  ion  acoustic  instabilities  in  infinite  homogeneous  plasmas.  The  problem  with  the  latter 
two  processes  is  that  they  are  slow  compared  to  the  laser  frequency  so  that  particle  simulations  get  very 
expensive.  For  instance,  the  simulation  of  Brillouin  backscatter  in  Ref.  11  ran  to  t  —  7  x  103/ fl 
(35,000  time  steps)  with  VJ  Vt  —  7.  However,  this  is  only  3.5  psec  for  a  Nd  laser  produced  plasma, 
and  the  total  length  corresponds  to  only  10  /*m,  meaning  a  density  gradient  scale  length  of  1-2  #xm. 
Clearly,  a  particle  simulation  for  more  realistic  times,  lengths  and  smaller  V0J  Ve  fields  would  be  astro¬ 
nomically  expensive.  Our  approach  sacrifices  the  detailed  description  of  the  instabilities.  Its  effect  is 
now  modeled  by  anomalous  transport  in  a  fluid  code.  However  we  gain  tremendously  in  the  range  of 


3 


LABORATORY  FOR  COMPUTATIONAL  PHYSICS  AND  PLASMA  PHYSICS  DIVISION 

parameters  that  can  be  studied.  For  instance,  for  a  70  psec  Nd  laser  pulse,  we  model  the  entire  under- 
dense  plasma  for  the  entire  pulse  duration  with  less  than  two  minutes  computing  time  with  a  TI-ASC. 

The  other  approach,  used  to  study  laser-plasma  interactions,  is  fluid  simulations  (characterized  by 
LASNEX  I2~IS).  Our  simulations  have  more  in  common  with  them  than  they  do  with  particule  simula¬ 
tions  but  there  are  still  important  differences.  Specifically  our  treatment  of  flux  limitation  is  quite 
different.  In  our  model,  a  large  thermal  flux  excites  an  instability.  This  instability  tends  to  limit  the 
flux,  and  gives  rise  to  several  other  anomalous  transport  effects.  All  of  these  transport  effects  can  be 
self-consistently  related  to  each  other  through  the  turbulence  level.  To  estimate  this  level  we  appeal  to 
previous  experiment,  simulation  and  theory.  Thus  flux  limits  are  derived  in  this  work,  not  specified. 
Also  this  work  incorporates  a  model  for  Brillouin  backscatter  and,  as  discussed  in  the  next  section,  our 
approach  to  electron  thermal  energy  transport  is  quite  different  from  that  used  in  LASNEX.15 

In  the  next  section  we  discuss  the  anomalous  physical  processes  in  our  simulation  models. 
Detailed  results  for  this  work  including  simulations  and  comparisons  with  experiment  may  be  found  in 
NRL  Memorandum  Report  4083. 
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Al.  Anomalous  Physical  Processes 

Ala.  Ion  Acoustic  Instability  Generated  by  Return  Current 

The  problem  of  the  ion  acoustic  instability  in  a  plasma  which  conducts  energy  has  been  discussed 
in  a  series  of  papers  which  deal  with  both  unmagnetized1, 16  and  magnetized2, 17  plasmas.  In  our  present 
work  we  confine  our  study  to  unmagnetized  plasmas.  The  basic  idea  is  that  if  the  electrons  conduct 
energy  but  carry  no  net  current,  the  flux  of  energetic  particles  in  the  direction  of  Q  (the  heat  flux)  must 
be  balanced  by  a  return  current  of  low  velocity  particles  in  the  opposite  direction.  This  return  current 
can  excite  ion  acoustic  waves  which  also  propagate  opposite  to  Q .  There  are  four  principal  physical 
effects  of  this  instability.  First  of  all  it  reduces  the  energy  flux,  or  in  other  words,  reduces  the  elec¬ 
tron  thermal  conduction.  Second,  it  gives  rise  to  an  anomalous  energy  exchange  between  electrons  and 
ions,  with  the  ions  gaining  energy  and  the  electrons  losing.  Third,  it  gives  rises  to  an  electric  field 
which  generates  the  return  current  (this  latter  process  is  not  directly  calculated  in  our  work  because  we 
deal  only  with  the  total  momentum,  not  with  the  momentum  of  each  species  separately).  These  three 
processes  are  discussed  in  detail  in  Ref.  16. 

Fourthly,  the  ion  density  fluctuations  give  rise  to  enhanced  absorption  of  the  laser  light.1,2,16,18 
Since  our  results  depend  upon  e<t>/Te  and  the  cone  angle  of  the  turbulent  spectrum  0,  it  is  important  to 
pick  these  parameters  as  accurately  as  possible.  Fortunately  we  can  draw  on  a  great  deal  of  experimen¬ 
tal  data19"26  and  particle  simulations27-33  on  current  driven  ion  acoustic  instabilities.  There  is  also  a 
great  deal  of  recent  theory34"38  which  examines  ion  trapping,  resonance  broadening  or  anomalous  tran¬ 
sport  and  which  supports  these  data.  Our  basic  assumption  is  that  an  ion  acoustic  instability  driven  by  a 
return  current  will  have  the  same  fluctuation  level  and  angular  spread  as  one  driven  by  a  real  current. 

Also  we  model  the  spectrum  with  a  single  k ,  A  choice  0  =  7t/3,  k  —  kD/ 2  and  —  0.1  seems  to  be 

‘e 

reasonably  consistent  with  the  data  which  we  have  cited. 

Let  us  finally  note  that  if  the  plasma  is  magnetized,  the  ion  wave  has  its  wave  vector  perpendicu¬ 
lar  to  both  B  and  V  Te,  so  the  wave  vector  can  be  parallel  to  the  electric  field  over  part  of  the  laser  spot. 
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Thus  not  only  can  the  absorption  be  enhanced,  but  also  there  is  no  need  to  make  any  assumption  con¬ 
cerning  angular  width  of  the  specturum.  Our  study,  which  neglects  B,  then  gives  rise  to  minimum 
anomalous  absorption. 

Alb.  Brillouin  Backscatter 

We  now  discuss  the  Brillouin  backscatter  instability  in  an  inhomogeneous  plasma.39  40,41  Since  the 
flow  velocity  in  the  underdense  plasma  is  generally  supersonic,  we  assume  that  the  group  velocity  of  the 
reflected  wave  and  acoustic  wave  are  in  the  same  (outward)  direction.  However,  even  if  the  group 
velocities  have  opposite  sign,  the  result  of  the  theory  is  not  different.39,42 

For  the  case  without  Landau  damping  of  the  sound  wave,  the  reflected  wave  travels  back  towards 
the  laser,  and  is  spatially  amplified  for  distances  less  than  a  length  Lc  based  on  the  the  growth  rate  of 
the  Brillouin  backscatter  instability.  But  if  the  sound  wave  is  strongly  damped  then  the  critical  length  is 
based  on  the  Landau  damping  rate  for  the  wave. 

If  this  linear  theory  is  used,  the  result  is  that  Brillouin  backscatter  is  an  extraordinarily  strong  pro¬ 
cess  and  it  dominates  the  dynamics  (that  is,  the  laser  light  is  almost  totally  reflected)  for  nearly  all  gra¬ 
dient  scale  lengths  and  incident  laser  light  intensities.  Particle  simulations  have  also  confirmed  that 
while  in  the  linear  regime,  Brillouin  backscatter  can  play  a  dominant  role  even  if  the  gradient  scale 
length  is  of  order  of  a  free  space  wavelength.11  Thus  a  crucial  problem  is  to  see  whether  any  process 
reduces  the  strength  of  Brillouin  backscatter  instability. 

We  consider  a  nonlinear  reduction  in  the  growth  rate  due  to  ion  trapping  by  the  ion  acoustic  wave. 
To  start,  we  will  show  that  as  the  ion  wave  grows,  it  is  forced  into  a  nonlinear  regime.  For  temporal 
growth,  without  Landau  damping,  each  photon  of  incident  wave  lost  produces  one  photon  of  reflected 
light  and  one  phonon  of  ion  acoustic  wave.  For  spatial  growth,  the  situation  is  different  because  the 
group  velocity  of  the  photon  is  much  greater  than  the  group  velocity  of  ion  acoustic  wave.42,44 


Now  let  us  consider  ion  trapping.  An  ion  will  trap  if  the  peak-to-peak  potential  energy  drop  of  the 
acoustic  wave  is  greater  than  the  difference  in  energy  between  the  wave  phase  velocity  and  V3  times 
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the  ion  thermal  energy  (assuming  a  waterbag  distribution  function  for  the  ions37).  Also,  for  trapping  to 
have  a  significant  effect  on  the  ion  distribution  function,  the  trapping  width  must  be  larger  than  the  ion 
thermal  speed.  Thus  we  have  a  double  condition  for  ion  trapping.  For  any  reasonable  parameters,  the 
saturation  amplitude  of  the  fluctuating  potential  of  the  ion  acoustic  wave  is  much  larger  than  that 
required  for  ion  trapping.  Therefore,  before  the  photon  can  be  strongly  backscattered,  the  ion  acoustic 
wave  is  driven  into  a  strongly  nonlinear  regime. 

To  describe  the  backscatter  in  this  regime,  we  invoke  recent  theoretical  work.37  45  Our  work  does 
not  utilize  an  enhanced  damping,  but  rather  a  modification  to  the  growth  rate  in  the  absence  of  damp¬ 
ing.  The  idea  is  that  there  is  some  power  input  to  the  ion  acoustic  wave.  In  the  linear  regime,  this 
power  input  causes  the  ion  acoustic  wave  to  grow  at  the  proper  linear  growth  rate.  However  once  the 
ions  begin  to  trap,  a  small  increase  in  amplitude  of  the  ion  acoustic  wave  traps  a  few  more  ions.  Since 
ions  abruptly  gain  a  tremendous  amount  of  energy  when  they  are  trapped,  the  power  input  goes  princi¬ 
pally  into  accelerating  ions  and  only  slightly  into  causing  wave  growth.  However  the  wave  does  not  stop 
growing  completely  because  the  power  input  is  still  there.  Nevertheless,  the  growth  rate  does  abruptly 
drop  once  ions  begin  to  trap.  The  calculation  of  the  growth  rate  is  given  in  Ref.  37  for  a  wave  which 
grows  in  time  in  an  homogeneous  plasma  and  in  Ref.  45  for  a  wave  which  grows  in  space  in  an  homo¬ 
geneous  but  bounded  plasma. 

In  summary,  our  non-linear  treatment  of  stimulated  Brillouin  backscatter  which  includes  pump 
depletion,  ion  heating  relies  on  reducing  the  linear  growth  rate  y0  and  the  coherent  scattering  length 
Lc.  It  is  significantly  different  from  other  SBS  treatments51  which  rely  mostly  on  increasing  Lc  until  it 
is  longer  than  the  plasma  size. 

Ale.  The  Problem  of  Energetic  Electrons 

Other  numerical  studies  of  this  problem,  principally  LASNEX  have  assumed  that  the  electron 
energy  is  transported  mostly  by  energetic  electrons  which  are  treated  separately.  The  basic  idea  is  that 
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the  classical  mean  free  path  is  much  longer  than  any  characteristic  scale  length.  For  instance  if  Te  - 
10  keV  and  n  -  1021  cm"3,  the  electron  mean  free  path  is  roughly  one  millimeter.  Thus  a  classical 
fluid  model  is  not  valid  for  these  electrons.  However  separate  transport  of  these  electrons  is  a  very 
complicated  procedure  and  to  our  knowledge,  there  is  no  generally  agreed  upon  procedure  for  it. 

The  approach  that  we  adopt  is  very  different.  We  use  a  single  fluid  model  and  justify  it  by  a  large 
anomalous  collision  frequency.  For  instance,  if  the  return  current  excites  an  ion  acoustic  instability 
with  fluctuation  amplitude  e<b/Te%  the  anomalous  collision  frequency  is  given  roughly  by  (e<fr/Te)2. 
Therefore,  if  e<t>/Te  —  10~l  the  electron  mean  free  path  for  a  10  keV  plasma  at  n  *  1021  cm-3  is 
roughly  2  microns.  Thus  a  single  fluid  approximation  is  now  much  more  reasoned-  Even  if  the 
absorption  process  gives  rise  to  nonthermal  electron  distributions  (for  instance  resonant  absorption 
almost  certainly  does8~u),  if  the  transport  of  these  energetic  electrons  is  inhibited,  a  single  fluid  theory 
is  probably  still  valid  everywhere  except  right  at  the  position  where  the  nonthermal  electrons  are 
accelerated.  At  these  points  the  thermal  temperature  we  calculate  will  be  some  average  of  the  superth- 
ermal  and  thermal  temperature  actually  produced. 

Thus  our  philosophy  is  to  explore  the  consequences  of  a  single  fluid  description  of  the  electrons 
and  see  where  it  leads.  For  one  thing  we  find  that  the  temperature  in  the  underdense  plasma  is  gen¬ 
erally  much  greater  than  the  thermal  electron  temperature  predicted  by  LASNEX.  Our  thermal  tem¬ 
perature  in  the  underdense  region,  in  fact  agrees  much  better  with  both  the  superthermal  electron 
energy  of  the  LASNEX  simulations  and  shown  by  experiment.  Of  course  there  will  still  be  runaway 
tails  on  the  electron  distribution  function,  but  now  they  perturb  rather  than  dominate  the  electron  tran¬ 
sport.  We  find  that  this  model  does  give  reasonable  results.  Another  great  advantage  is  that  it  is  a  very 
straightforward  theory  and  is  very  economical  to  study. 

A2.  Conclusions 

We  conclude  that  fluid  simulations  with  anomalous  transport,  is  an  accurate  and  cost  effective  way 
to  describe  the  absorption,  backscatter  and  flux  limitation  in  a  laser  produced  plasma.  Our  calculations 
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give  quantitatively  correct  absorption  of  laser  light  over  five  orders  of  magnitude  in  irradiance,  for  short 

pulses,  long  pulses,  and  double  structured  pulses.  It  has  been  shown  that  reducing  the  laser  wavelength 
is  favorable  to  absorption  of  laser  light  by  targets  (although  it  might  not  be  overall  more  efficient). 
Another  important  result  uncovered  by  this  work  is  that  the  physics  of  absorption  is  different  for  short 
and  long  pulses  at  high  intensities  (d»k2  >  1014  W/c m2~n2).  This  has  been  confirmed  recently  in 
experiments.56 
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Bl.  Energy  Conversion  of  Light  Ion  Beams  by  Ablative  Acceleration  of  Thin  Foils 

Light  ion  beams  provide  a  power  source  of  longer  pulse  duration  than  that  typical  of  laser  power 
sources.  This  source  in  turn  calls  for  larger  size,  thin  shell  pellet  designs  for  inertial  confinement 
fusion.  The  optimum  thickness  of  the  foil  for  conversion  of  beam  energy  to  kinetic  energy  of  the  tar¬ 
get  payload  depends  on  making  the  foil  sufficiently  thick  to  stop  the  beam  completely  and  maintain  tar¬ 
get  integrity  during  beam  absorption  while  keeping  the  mass  of  the  payload  to  a  minimum,  maximizing 
the  final  velocity  achieved. 

We  report  here  the  results  of  a  numerical  investigation  of  the  accelerations  of  a  variety  of 
thicknesses  of  thin  gold  foils  by  a  1.4  Mev  proton  beam.  The  proton  beam  is  simulated  by  a  time- 
dependent  monoenergetic  beam  of  parallel  protons.  The  time  history  of  the  energy  and  current  per 
square  centimeter  of  the  beam  shown  in  Fig.  Bl.l  corresponds  to  the  beam  produced  experimentally  in 
the  Gamble  II  facility  at  NRL.  It  is  modeled  by  linear  segments  with  a  maximum  beam  energy  of  1.4 
Mev  and  maximum  beam  current  of  450  kA/cm2  this  corresponds  to  a  maximum  beam  power  of  6.3 
x  10* 1  watts/cm2.  The  current  pulse  length  is  65  ns  which  results  in  a  total  energy  in  the  beam  over 
the  pulse  of  19.4  KJ/cm2. 

The  target  consists  of  a  uniform  plane  gold  foil  with  a  parallel  proton  beam  impinging  from  the 
right.  The  convention  will  be  used  that  the  front  side  of  the  foil  is  the  side  on  which  the  proton  beam 
impinges.  The  numerical  simulations  are  done  for  a  uniform  slab  of  thickness  ranging  from  7.5/x  to 
2Q u-  The  1.4  Mev  proton  beam  has  a  range  of  I<V  in  the  cold  gold  foil  so  that  the  foil  thickness 
ranges  from  subrange  to  twice  the  beam  range  in  the  foil. 

The  numerical  simulations  were  carried  out  using  a  quasi- Lagrangian  one-dimensional  hydro  code 
in  slab  geometry.  The  essential  features  of  the  code  are  outlined  in  Fig.  B1.2.  The  code  solves  the 


14 


NRL  MEMORANDUM  REPORT  *369 

conservation  equations  for  mass,  momentum  and  energy  with  a  flux  conserving  numerical  scheme. 
The  interfaces  with  the  free  boundaries  on  either  side  are  treated  in  a  fully  Lagrangian  manner.  The 
regions  between  the  two  outer  interfaces  are  treated  using  the  moving  zone  flux-corrected  transport 
(FCT)  algorithms.  The  zones  are  moved  in  such  a  way  that  the  mass  per  unit  zone  is  approximately 
constant.  The  code  incorporates  a  single  fluid,  two-temperature  (electron  and  ion)  model  with  an 
independent  electron  energy  transport  equation.  The  ion  temperature  is  found  by  subtracting  the  elec¬ 
tron  and  kinetic  energies  from  the  total  energy.  Classical  electron  and  ion  thermal  conductivities  are 
treated  in  an  implicit  conduction  routine  along  with  the  equilibration  terms.  In  the  runs  to  be  described 
here  the  electrons  and  ions  are  constrained  to  have  the  same  temperature  assuming  local  thermo¬ 
dynamic  equilibrium.  Radiation  transport  is  treated  using  the  radiation  conduction  approximation  in  the 
optically  thick  region  where  local  thermodynamic  equilibrium  is  assumed.  The  optically  thin  regions  are 
treated  using  a  volume  of  emission  coefficient.  The  transition  from  thick  to  thin  is  assumed  to  take 
place  one  optical  depth  from  each  of  the  free  surfaces.  The  brightness  temperature  is  calculated  based 
on  the  integrated  radiation  from  each  side.  The  last  thick  cell  on  each  side  is  limited  by  black  body 
radiations  from  that  cell. 

The  equation  of  state  of  the  hot  gold  plasma  is  calculated  from  curve  fits  to  the  Saha  equilibrium 
calculations  for  gold.  At  low  temperatures  the  solid  material  properties  of  the  gold  are  included  from 
empirical  curve  fits  to  the  cold  compression  curves  and  using  a  Griineisen  coefficient  for  the  thermal 
component.  The  equation  of  state  is  matched  at  the  low-temperature,  high  density  limit  to  the 
Thomas-Fermi-degenerate  electron  limit  and  at  high  temperatures  and  low  densities  to  the  Saha  equili¬ 
brium  limit. 

The  ion  beam  deposition  is  treated  in  two  parts.  The  slowing  down  due  to  bound  electrons  is 
treated  using  an  empirical  experimentally  determined  range  energy-relation  to  determine  the  slowing 
down  rate  as  a  function  of  beam  energy  and  target  density.  This  determines  the  cold  slowing  down  for¬ 
mulae.  To  this  is  added  a  correction  as  the  gold  becomes  ionized  for  the  bound  ionization  potential  and 
due  to  the  slowing  down  due  to  free  electrons  which  are  treated  classically. 
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The  results  may  be  summarized  in  Figs.  B1.3  and  B1.4.  In  Fig.  B1.3  we  show  the  asymptotic  dis¬ 
tribution  of  energy  of  the  gold  foil  as  a  function  of  the  foil  thickness  (after  it  has  expanded  to  low  pres¬ 
sure).  This  is  the  measure  of  how  the  beam  energy  will  be  distributed  after  the  foil  has  accelerated  to 
maximum  velocity.  The  payload  kinetic  energy  is  defined  as  the  energy  of  all  the  gold  material  moving 
in  the  direction  away  from  the  beam  source.  It  shows  a  broad  maximum  in  the  range  of  lO-15/i  gold 
thickness.  Below  10/i  not  ail  of  the  beam  energy  is  absorbed  as  it  is  insufficiently  thick  to  stop  the 
beam.  Above  15/a  it  begins  to  fall  off  as  the  increased  mass  of  the  payload  reduces  the  final  maximum 
velocity  obtained.  The  total  momentum  delivered  to  the  payload  remains  nearly  constant  as  a  function 
of  foil  thickness  once  the  beam  is  totally  absorbed,  thus  the  increased  mass  results  in  a  lower  payload 
velocity.  The  internal  energy  includes  both  thermal  and  ionization  energy  and  increases  only  slowly 
with  increased  foil  thickness  representing  about  12%  of  the  total  energy.  The  radiated  energy  increases 
at  first  as  more  of  the  beam  energy  is  absorbed  but  then  decreases  when  the  foil  becomes  thicker  than 
the  beam  range  since  it  radiates  from  one  side  only  until  the  thermal  wave  bums  through  to  the  back 
surface.  This  takes  longer  with  the  thicker  foil  and  after  the  beam  shuts  off  the  foil  cools  due  to  adia¬ 
batic  expansion  before  burning  through,  thus  radiating  less.  Most  of  this  energy  is  recovered  in  kinetic 
energy  of  the  foil  material  as  can  be  seen  by  the  rapid  rise  in  total  kinetic  energy.  However,  most  of 
the  kinetic  energy  is  carried  off  with  the  ablating  foil  material,  hence  the  broad  flat  peak  in  the  payload 
kinetic  energy.  The  increased  energy  available  from  the  lowered  radiated  energy  balances  the  decrease 
in  payload  kinetic  energy  due  to  the  increasing  mass  of  the  payload.  The  result  is  a  broad  peak  where 
a  broad  peak  where  about  27%  of  the  beam  energy  has  been  converted  to  payload  kinetic  energy  mov¬ 
ing  at  an  average  velocity  of  about  2  x  10*  cm/sec. 

In  Fig.  B1.4  we  show  the  brightness  temperature  of  the  front  and  back  sides  of  the  foil.  Also 
shown  is  the  peak  temperature  anywhere  in  the  foil  as  a  function  of  foil  thickness.  This  temperature 
rises  slowly  with  increased  foil  thickness  as  the  result  of  the  more  massive  foil  remaining  intact  longer 
reducing  the  cooking  from  adiabatic  expansion.  The  brightness  temperature  is  a  measure  of  the  total 
radiation  from  an  optically  thick  plasma.  This  figure  demonstrates  the  dramatic  drop  in  radiation  from 
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the  back  surface  as  the  foil  thickness  exceeds  the  range  of  the  ion  beam.  When  the  ion  beam  no  longer 
penetrates  to  the  back  surface  no  radiation  takes  place  until  a  thermal  wave  from  the  hot  plasma  inte¬ 
rior  burns  through  to  the  surface.  As  the  foil  becomes  thicker  this  takes  a  longer  time  and  the  plasma 
interior  begins  to  cool  by  adiabatic  expansion  and  radiation  on  the  from  side  before  the  wave  can  burn 
through,  thus  the  back  surface  will  radiate  at  a  much  lower  temperature  and  for  a  shorter  period  of 
time. 


In  Figs.  B1.5  to  B1.7  we  show  the  velocities  and  temperatures  of  the  front  and  back  surface  as  a 
function  of  time.  Plotted  is  the  velocity  of  the  leading  edge  of  :he  front  side  and  the  negative  of  the 
back  side  velocity.  Also  shown  is  the  mean  velocity  of  the  payload  which  is  defined  as  the  total 
momentum  of  material  moving  in  the  direction  away  from  the  beam  divided  by  the  mass  of  this 
material.  The  closer  this  velocity  is  to  the  back  side  edge  the  more  nearly  the  payload  is  being 
accelerated  as  an  intact  foil.  In  the  7.5m  case  shown  in  Fig.  B1.5  the  ion  beam  penetrates  to  the  back 
side  edge  and  accelerates  it  to  high  velocity.  If  we  look  at  the  temperature  of  the  back  side,  it  does  not 
begin  its  rapid  rise  until  about  1 2  ns.  This  is  the  result  of  using  the  ramp  profile  in  time  for  the  beam 
energy.  The  early  portion  of  the  beam  does  not  have  sufficient  energy  to  penetrate  the  foil.  Note  that 
the  peak  temperature  occurs  at  about  45  ns  which  is  the  time  of  peak  ion  beam  power.  In  Fig.  B1.6  we 
see  the  results  for  a  foil  thickness  of  i 2.5m  which  is  somewhat  greater  than  the  proton  beam  range. 
Note  that  for  the  first  20  ns  the  leading  edge  and  average  velocities  are  nearly  the  same  since  the 
accelerated  leading  edge  remains  a  solid.  At  about  25  ns  a  thermal  wave  from  the  hot  interior  reaches 
the  back  surface  accompanied  by  a  rapid  increase  in  the  brightness  temperature  and  a  rapid  acceleration 
of  the  back  edge  as  the  surface  vaporizes  and  blows  off.  This  rapid  expansion  results  in  a  cooling  which 
can  be  seen  by  the  lowering  of  the  radiation  temperature.  At  later  times  a  second  thermal  wave  catches 
up  with  the  leading  edge  causing  another  rise  in  the  surface  temperature.  In  addition,  as  the  density  of 
the  leading  edge  drops,  the  outer  plasma  becomes  optically  thin  and  radiation  from  the  hotter  interior 
escapes.  Note  also,  as  in  the  7.5m  case,  most  of  the  acceleration  takes  place  while  the  ion  beam  is  irra¬ 
diating  the  target.  Since  the  foil  becomes  completely  vaporized  the  leading  and  trailing  edges  both  have 
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much  greater  velocity  than  the  mean  payload  velocity.  For  the  20m  foil  case  shown  in  Fig.  B1.7  the  foil 
is  sufficiently  thick  that  the  thermal  wave  does  not  arrive  at  the  front  surface  until  after  the  proton 
beam  has  ended  and  is  at  a  much  lower  temperature.  Note  that  the  mean  velocity  and  the  leading  edge 
velocity  remain  nearly  the  same  indicating  acceleration  at  near  solid  density.  A  weak  shock  arrives  at 
the  leading  edge  at  about  30  ns  but  is  insufficient  to  vaporize  the  gold  material  yielding  only  an 
increased  acceleration  of  the  leading  edge.  The  thermal  wave  which  finally  reaches  the  leading  edge  at 
about  75  ns  is  much  colder  than  the  others  since  the  plasma  has  begun  to  cool  before  it  reaches  the 
leading  edge.  The  front  side  temperature  and  leading  edge  velocity  are  relatively  the  same  for  each 
case.  This  yields  the  result  that  the  momentum  transferred  to  the  payload  is  nearly  constant.  The  velo¬ 
city  and  kinetic  energy  decrease  as  the  mass  of  the  payload  increases. 

The  results  of  these  calculations  indicate  that  at  least  at  these  power  levels  an  ion  foil  can  be 
accelerated  with  greater  than  25%  efficiency  in  conversion  of  beam  energy  to  kinetic  energy  of  the  mov¬ 
ing  foil.  The  acceleration  is  ablative  as  can  be  seen  by  the  fact  that  most  of  the  foil  acceleration  takes 
place  while  the  beam  power  is  on,  and  the  results  indicate  that  the  rocket  model  where  the  momentum 
transfer  is  proportionate  to  the  ablated  mass  and  ablation  velocity  is  applicable  for  the  ion  beam 
acceleration.  Since  the  ion  beam  penetrates  approximately  a  constant  mass  of  material  for  a  uniform 
beam  energy,  the  results  indicated  that  most  of  the  ablation  mass  is  contained  in  the  region  the  beam 
has  penetrated  and  little  additional  mass  is  ablated  from  the  penetration  of  the  thermal  wave. 
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Fig.  B1 . 1  —  Ion  Beam  current  and  voltage  profiles  as  a  function  of  time.  Zero  lime  corresponds  to  start  of  beam  deposition. 
Diagram  at  right  shows  schematic  of  beam  deposition  profile  in  cold  gold  foil. 
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FCT  ONE-DIMENSIONAL  HYDROCODE 


•  Solves  fully  conservative  equations  with  a  flux- 
conserving  numerical  scheme  in  either  slab,  cylin¬ 
drical,  or  spherical  geometry. 

•  Interfaces  between  materials  and  free  boundaries 
are  fully  Lagrangian  with  arbitrary  numbers  of 
interfaces. 

•  Regions  between  interfaces  are  treated  with  moving 
zone  flux-corrected  transport  (FCT). 

•  Two- temperature  model  (Te,  Tj)  is  used  with  implicit 
thermal  conduction  and  equilibration  (classical  coeffi¬ 
cients  with  flux  limiting). 

•  Radiation  transport  uses  conduction  approximation  in 
optically  thick  regions  and  volume  emission  in  thin 
regions  with  transition  at  one  optical  depth  from 
surface. 

•  Equation-of-state  treats  solid  material  properties 
and  Ferm-degenerate  electrons  at  high  densities  with 
ionization  state  determined  by  Saha  equilibrium. 

•  Ion  beam  deposition  uses  empirical  range  energy 
relation  for  bound  electrons  and  classical  slowing 
down  for  free  electrons. 

Fig.  B1.2  —  Essential  features  of  the  Quasi-Lagrangian  one-dimensional  hydro  code  in  slab  geometry. 
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ION-FOIL  ENERGY  BUDGET 


FOIL  THICKNESS  (p) 

Fig.  Bt  .3  —  Distribution  of  energy  after  foil  has  expanded  for  150  ns.  Payload  kinetic  energy  peaks  at  a  foil  thickness 
that  approximately  corresponds  to  the  ion  beam  range  in  cold  gold. 
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ION-FOIL  TEMPERATURE 
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Fig.  B1.5  —  Velocity  and  temperature  of  front  and  rear  surfaces.  V corresponds  to  the  average  velocity  of  the  payload  as 

distinct  from  the  back  surface  edge  velocity. 
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by  expansion  of  buck  surface  after  the  thermal  wave  reaches  it 
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TIME  (ns) 

Fig.  BI.7  —  Velocity  and  temperature  of  front  and  rear  surfaces  The  average  velocity  and  back  side  velocity  are  nearly  the  same 
as  long  as  the  foil  remains  solid  until  it  is  heated  sufficiently  to  vaporize. 
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B2a.  Spherical  Shock  Implosion  Dynamics  and  Stability 

One  of  the  critical  limitations  to  achieving  a  high  degree  of  compression  in  a  spherical  implosion 
is  the  degree  of  symmetry  that  can  be  maintained  during  the  implosion  process.  Several  procedures  are 
available  to  help  analyze  the  issues  of  the  stability  and  symmetry  in  implosions  involving  both  analytic 
and  numerical  treatment.  The  simplest  of  these  involves  linearized  perturbation  analysis  of  some  zero 
order  spherically  symmetric  solution  to  determine  if  the  solution  is  stable.  If  it  is  not,  multidimensional 
calculations  are  required  to  estimate  how  large  the  perturbations  become  and  if  there  is  a  saturation 
mechanism.  Rayleigh-Taylor  modes  in  a  laser-driven  ablation  layer  have  been  analyzed  extensively  at 
NRL  in  this  manner. 

The  zero  order  solution  may  be  found  either  analytically  using  some  self-similar  solution,  or 
numerically.  The  perturbations  also  may  be  investigated  either  analytically  or  numerically  leading  to  the 
so-called  "piggy-back"  methods.  A  complementary  approach  to  finding  the  growth  of  linearized  unstable 
modes  is  the  direct  numerical  solution  of  the  multi-dimensional  fluid  equations.  In  this  approach  the 
solutions  are  found  without  perturbation  analysis  at  the  expense  of  a  rather  large  hydrodynamics  calcu¬ 
lation  to  be  performed  for  each  set  of  parameters  to  be  investigated.  Modes  of  a  few  zone  sizes  are  the 
smallest  that  can  be  investigated  in  a  detailed  multidimensional  simulation.  This  approach  requires 
rather  fine  gridding  and  expensive  calculations  to  look  at  short  wavelength  phenomena. 

One  of  the  important  issues  in  understanding  imploding  systems  is  the  stability  of  a  converging 
shock  wave.  This  shock  might  be  used,  for  instance,  not  only  to  compress  the  fuel,  but  also  to  provide 
the  required  heating  to  create  a  central  ignition  region.  The  final  temperature  achievable  will  depend  on 
how  spherical  the  shock  wave  remains  during  the  collapse  process  and  the  shape  of  the  shock  at  the 
time  of  reflection. 

The  motion  of  a  converging  shock  wave  can  be  computed  with  great  accuracy  by  considering  only 
the  changes  in  the  physical  variables  across  the  shock  front,  ignoring  the  motion  behind  the  shock 
front.  The  shock  front  moves  normal  to  the  shock  surface  so  it  may  be  treated  as  locally  one- 
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dimensional  flow  down  a  channel  whose  boundaries  are  determined  by  the  trajectories  of  the  shock 
front.  These  trajectories  form  imaginary  ray  tubes  whose  cross-sectional  area  may  be  related  to  the 
Mach  number  by  an  equation  developed  by  Chester-Chisnell-Whitham  (CCW  approximation)1-3 
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This  approximation  allows  a  great  deal  of  savings  in  computational  time  and  storage  since  one  needs 
only  store  the  location  and  Mach  number  of  the  shock  front  and  compute  the  trajectory  of  the  shock 
front  rather  than  the  entire  flow  field.  The  problem  can  be  reduced  to  the  calculation  of  a  set  of  ordi¬ 
nary  differential  equations  for  the  trajectories  of  a  finite  number  of  points  located  along  the  shock  front 
and  for  their  associated  time-dependent  Mach  numbers. 


A  code  based  on  this  approximation  has  been  written  and  enables  us  to  investigate  a  wide  range  of 
stability  and  symmetry  issues  related  to  shock  collapse.  For  instance,  there  is  the  question  of  how  large 
the  initial  perturbations  may  be  even  for  stable  modes  since  the  perturbation  amplitude  may  not  go  to 
zero  as  fast  as  the  average  radius.  The  concept  of  stability  in  the  sense  of  a  growing  or  decaying  pertur¬ 
bation  does  not  make  much  sense  for  this  particular  problem.  Given  a  set  of  wavelengths,  Mach 
number,  and  initial  radius,  can  the  collapse  of  the  shock  be  timed  so  that  the  mode  amplitude  of  an 
oscillating  perturbation  is  zero  at  the  time  of  collapse?  Is  there  a  cutoff  mode  number  below  which  the 
wave  front  becomes  more  rather  than  less  spherical  during  collapse? 

In  order  to  address  these  questions  we  have  developed  a  code  similar  to  that  of  Fong  and 

Ahlbom4.  The  motion  of  the  shock  front  is  described  by  K  mesh  points  whose  equation  of  motion  is 

given  by: 

xrl  -  xa + Mr*  cos  -  n  «a> 

-  y™  +  A sin  {Ng+^ait**1  -  tn)  (2b) 

where  (X£f  YU)  is  the  Cartesian  coordinate  of  the  fcth  mesh  point  at  the  nth  time  interval;  Mg**  is  the 

Mach  number  at  N£*'h  is  the  angle  between  the  shock  and  the  axis,  a  is  the  sound  speed 

ahead  of  the  shock,  and  t  the  time.  The  n+'h  values  are  formed  by  advancing  the  equations  for  half  a 
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timestep  to  make  the  overall  algorithm  second  order.  The  angle  factors  are  evaluated  by  averaging  the 


angles  formed  by  joining  the  adjacent  mesh  points  with  straight  lines. 
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The  mach  number  is  given  from  the  CCW  approximations  as: 
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where  the  area  of  the  ray  tube  between  mesh  points  is  given  by 
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This  gives  a  second-order  space  and  time  centered  algorithm  which  is  stable  subject  to  a  timestep  limit 
analogous  to  the  Courant  condition  for  the  one-dimensional  fluid  equations.  That  is 

M  <  b  mm  (^fA’  r"k--  (7) 

MU  a 

where  b  is  a  numerical  factor  less  than  one  chosen  to  prevent  nonlinear  instabilities  from  developing. 
We  have  used  values  of  b  -  .05.  This  algorithm  allows  the  propagation  of  shock-shocks  (discontinui¬ 
ties  in  the  slope  and  mach  number  along  the  shock  front  predicted  by  the  Whitham  theory)  in  either 
direction  since  the  scheme  is  centered  and  symmetric.  It  is  not  clear  that  the  algorithm  proposed  by 
Fong  and  Ahlborn  accomplishes  this.  As  in  the  Fong  and  Ahlborn  procedure  we  found  it  necessary  to 
redistribute  the  grid  points.  This  was  done  not  for  stablity  but  to  prevent  unacceptably  short  timesteps 
as  the  mesh  points  crowd  together  near  the  shock-shock  regions.  The  algorithm  used  here  diffuses  the 
grid  points  a  small  amount  parallel  to  the  shock  front  in  such  a  way  as  to  make  the  distances  between 
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grid  points  more  nearly  equal.  This  plays  a  role  similar  to  that  of  artificial  viscosity  in  the  one- 
dimensional  fluid  algorithm  or  surface  tension  at  an  interface. 


As  a  test  of  the  accuracy  of  the  numerical  procedure  described  above,  we  compared  the  results  of 
the  calculation  with  that  predicted  by  the  self-similar  solutions  for  the  collapse  of  an  infinitely  strong 
shock  due  to  Guderley.5  The  self-similar  solution  predicts  a  shock  position  given  by  R  *  A  (  —  /)“ 
with  with  a  —  .717  for  a  y  —  1.4  gas  given  by  Guderley.  This  gives  a  shock  Mach  number  as  a  func- 

0-1  «-l 

tion  of  radius  R  **  —  aA  (  —  t)a~]  —  R  a  .  Thus  M  —  R  a  —  /?“  3947.  The  CCW  approximation 
for  large  M  gives 
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This  gives  M  —  R  **  /?“  3941.  This  agrees  with  the  analytic  Guderley  solution  to  within 

.15%.  In  Fig.  B2a.l  we  show  the  results  of  these  models  for  the  spherical  convergence  of  a  shock  with 
an  initial  mach  number  of  57  at  a  radius  of  .1 12  cm.  the  solid  line  is  the  analytic  result  of  the  Guderley 
solution.  The  triangles  are  the  numerical  results  of  the  CCW  approximations  using  the  code  described 
above  with  25  mesh  points  around  a  circle  in  cylindrical  coordinates  so  that  the  result  describes  a  spher¬ 
ical  implosion.  Comparing  the  results  along  the  axis  to  those  perpendicular  to  the  axis,  the  shock 
remains  spherical  to  better  than  1%  during  the  entire  implosion  and  the  mach  number  agrees  with  the 
Guderley  solution  with  less  than  .5%  error.  Since  the  CCW  approximation  breaks  down  as  the  shock 
reflects  from  the  origin,  it  is  important  to  calibrate  as  well  a  method  capable  of  computing  past  the 
implosion  into  the  explosion  phase.  The  vertical  bars  are  the  results  computed  using  the  one¬ 
dimensional  hydrocode  FCT1D.  The  length  of  the  vertical  bar  represents  the  uncertainty  in  the  results 
of  the  pressure  behind  the  shock  due  to  an  uncertainty  in  its  location  on  a  finite  grid  and  the  steep  pres¬ 
sure  gradient  behind  the  shoe*.  These  difficulties  are  one  of  resolution  and  can  be  refined  by  increased 
resolution  at  the  expense  of  computer  time  and  storage. 


As  a  further  check  on  the  accuracy  of  the  FCT1D  solution,  the  ratio  of  the  pressure  at  a  given 
location  to  the  pressure  behind  the  first  incident  shock  at  the  location  is  given  as  a  function  of  the  ratio 
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of  time  to  the  time  of  collapse  of  that  shock  and  compared  to  the  analytic  solution.  Since  the  problem 
is  self-similar,  this  plot  should  be  the  same  for  all  physical  locations  in  space.  In  Fig.  B2 a.2  we  plot 
P/Pt  vs.  (t  -  tc)Kt0  —  tc)  where  the  time  of  collapse  tc  has  veen  taken  to  be  zero.  The  solid  line 
represents  the  analytic  solution  given  by  Guderiey.  The  horizontal  dashed  line  represents  the  theoreti¬ 
cal  jump  of  26  behind  the  reflected  shock.  The  vertical  bars  with  the  x’s  at  the  end  represent  the 
spread  in  the  values  found  from  the  numerical  results.  The  uncertainty  in  the  pressure  of  the  numeri¬ 
cal  results  is  due  to  both  an  uncertainty  in  the  exact  location  of  the  incoming  shock  used  to  determine 
the  reference  shock  location  (due  to  the  finite  width  of  the  shock)  and  an  uncertainty  in  the  exact  pres¬ 
sure  behind  this  shock  (due  to  the  steep  gradient  in  pressure  behind  the  shock  and  the  finite  shock 
width).  The  results  follow  the  analytic  solution  within  the  accuracy  with  which  they  can  be  interpreted 
from  the  numerical  solution.  This  gives  us  some  confidence  that  the  numerical  procedure  that  we  are 
using  can  accurately  determine  the  solution  for  the  shock  collapse  problem  and  subsequent  reexpansion 
with  no  adjustable  parameters  needed. 

Using  the  C-C-W  code  we  investigate  the  stability  of  a  perturbed  shock  wave.  Figures  B2a.J  and 
B2a.4  summarize  pictorially  the  temporal  evolution  of  an  8th  order  Legendre  polynomial  perturbation 
to  a  spherical  shock  with  amplitude  of  5%  and  10%,  respectively.  As  the  shock  collapses  toward  the  ori¬ 
gin,  the  amplitude  of  the  10%  perturbation  seems  to  grow  relative  to  the  radius.  After  implosion  by  a 
factor  of  5:1,  the  10%  perturbation  goes  unstable  whereas  the  5%  perturbation  is  still  behaving  within 
acceptable  bounds.  Shown  is  the  shape  of  the  shock  wave  as  it  collapses  from  right  to  left  with  increas¬ 
ing  time. 

We  also  have  investigated  how  symmetric  the  driver  must  be  in  order  to  produce  an  imploding 
shock  which  is  nearly  spherical.  Using  the  Chester-Chisnell-Whitham  (CCW)  approximation  we  have 
studied  a  wide  range  of  stability  and  symmetry  issues  related  to  shock  collapse.  These  yield  a  set  of 
wavelengths,  mach  numbers  and  initial  ratio  for  which  the  collapse  of  the  shock  can  be  timed  so  that 
the  amplitude  of  any  particular  oscillating  perturbation  is  near  zero  at  the  instant  of  collapse.  Superpo¬ 
sition  of  (individually)  unstable  modes  can  be  stabilized  nonlinearly.  We  have  also  determined  as  a 
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function  of  mode  number  and  initial  perturbation  amplitude  the  point  at  which  the  growth  becomes 
nonlinear.  The  stability  has  been  investigated  in  a  special  case  of  the  Guderley  solution  for  self-similar 
spherical  shocks  imploding  into  a  medium  having  a  power-law  density  profile.  The  linearized  equations 
of  motion  yield  an  exact  solution  which  predicts  power-law  growth.  This  result  agrees  with  that  found 
from  the  CCW  approximation,  which  has  been  extended  to  nonuniform  densities. 
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Spherical  Shock  Collapse 


Radius 


Fig.  B2a  I  —  Comparison  of  shock  mach  number  vs  radius  for  three  models  <0  Self*similar  Guderley  solution.  (2)  Numerical 
solution  of  the  C  C-W  approximation  equations  (31  Numerical  solution  of  I  D  fluid  equations  with  FCTID 
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SPHERICAL  SHOCK  COLLAPSE 


=  TIME  RELATIVE  TO  SHOCK  COLLAPSE 
l*®1  TIME,  t0 

Fig  B2a  2  —  Comparison  of  analytic  solution  due  to  Gudertey  with  the  results  of  the  FCT1D  code  for  the  pressure  time  history 

of  a  fixed  location  during  the  shock  collapse  and  reexpansion 


PERTURBED  SPHERICAL  SHOCK  COLLAPSE 


Fig.  B2a.3  —  Solution  using  the  C-C  W  approximation  lor  shock  collapse  of  a  spherical  shock  with  an  8th  order  Legendre  polyno¬ 
mial  perturbation  of  5%  amplitude  The  perturbation  becomes  smaller  at  the  same  rale  as  the  average  radius  remaining  within 
acceptable  bounds. 


PERTURBED  SPHERICAL  SHOCK  COLLAPSE 


Fig.  B2a  4  —  Solution  using  the  C-C-W  approximation  for  shock  collapse  with  a  10%,  8th  order  Legendre  polynomial 
perturbation  to  a  spherical  shock  The  perturbation  seems  to  grow  relative  to  the  average  radius  indicating  an  unstable 
behavior 
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B2.b  Modeling  of  the  Laser  Ablation  of  Thin  Foils 

It  has  been  demonstrated  in  a  series  of  experiments  at  NRL1  and  elsewhere2  that  the  use  of 
longer  duration  and  lower  intensity  laser  pulses  is  useful  for  achieving  successful  ablative  accelerations 
in  laser-pellet  fusion  studies.  At  NRL,  a  modified  version  of  the  Pharos  II  laser  system3  was  used  to 
irradiate  thin  foil  targets  of  plastic  and  aluminum.  Diagnostics  such  as  interferometry,  shadowgraphy, 
Doppler  sounding,  calorimetry,  and  charge  collection  were  used  to  deduce  both  target  and  ablation 
velocity,  mass  loss,  hydrodynamic  efficiency  and  ablation  pressure.1  These  NRL  measurements  provide 
a  reasonably  accurate  and  complete  macroscopic  picture  of  the  ablation  process,  but  quantitative  micros¬ 
copic  understanding  cannot  be  inferred. 

We  have  analyzed  the  NRL  experiments  numerically  in  order  to  quantify  the  controlling  physical 
principles  at  work  in  this  ablative  acceleration  and  to  study  the  sensitivity  of  the  macroscopic  results  to 
plausible  variations  in  the  initial  transients  and  the  electron  thermal  conduction  transport.  The  ablative 
acceleration  that  results  when  the  longer,  lower  intensity  pulses  are  used  has  been  the  subject  of  con¬ 
tinuing  concern.4 

The  ID  computer  model  used  combines  a  new  hydrodynamic  package  (ADINC)5  with  a  classical 
plasma  thermal  conduction  routine  (DIFFU1).  The  computational  grid  is  Lagrangian,  variably-spaced, 
and  rezonable.  Since  the  grid  is  Lagrangian  and  the  actual  problem  we  wish  to  consider  is  time  depen¬ 
dent,  the  initial  cell  sizes  can  change  in  time.  In  order  to  prevent  cells  from  becoming  either  too  large 
or  too  small,  a  rezoning  package  has  been  included  in  our  program.  It  is  capable  of  combining  neigh¬ 
boring  cells  while  conserving  both  energy  and  mass.  It  also  can  insert  a  new  cell  boundary  in  regions 
where  either  the  cells  are  becoming  too  large  or  temperature  gradients  are  to  steep.  This  rezoning 
allows  us  to  use  relatively  large  time  steps  (on  the  average  about  8  psec)  although  the  flow  velocities 
are  high  in  the  blow-off  region  (about  6  x  107  cm/sec)  and  the  cell  sizes  small  in  the  ablation  region 
(about  0.1  /a). 
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The  2D  code  is  FAST2D.6  It  uses  a  sliding  Eulerian  grid,  the  Flux-Corrected  Transport  (FCT) 
algorithm  for  solving  the  hydrodynamic  equations,  and  a  two-dimensional  classical  plasma  thermal  con¬ 
duction  routine. 

For  both  codes,  the  thermal  conduction  coefficient  has  the  classical  nonlinear  Pn  temperature 
dependence.  We  use  the  standard  Braginskii  coefficient7  and  assume  ail  plasmas  are  fully  ionized. 
Also,  for  the  ID  and  2D  calculations,  the  laser  flux  is  introduced  as  an  energy  source  term  in  the  heat 
equation  and  is  deposited  into  the  two  cells  about  the  critical  density.  The  grid  is  initially  uniform  in 
the  region  of  the  foil.  Inside  and  outside  of  this  region  the  grid  is  stretched.  On  the  inside  (opposite 
from  the  laser)  this  stretched  region  is  filled  with  a  very  low  density  plasma  having  the  same  adiabat  as 
the  target  plasma.  On  the  outside,  beyond  the  critical  density,  both  the  density  and  pressure  drop  off 
exponentially. 

The  initial  density,  pressure,  and  temperature  profiles  are  generated  from  the  analytical,  quasi¬ 
static  equilibrium  solutions  of  J.  H.  Orens.8  We  have  shown  that  these  analytic  solutions  provide  an 
adequate  steady  state  within  our  ID  computer  code.  This  is  done  by  initializing  test  runs  with  these 
profiles  and  then  proceeding  for  i,000  timesteps  at  constant  laser  intensity.  At  the  end  of  these  tests 
only  relatively  minor  changes  in  the  initial  profiles  have  occurred. 

Results  of  our  calculations  for  the  standard  experimental  case  of  a  15 fx  thick  plastic  (CH)  foil  irra¬ 
diated  by  a  laser  pulse  with  approximately  a  3  nsec  fwhm  and  a  peak  intensity  of  1013  w/cm2  are  sum¬ 
marized  in  Table  B2b.l,  Here,  A/0  is  the  initial  mass  of  the  foil,  is  the  mass  loss,  r\h  is  the  hydro- 
dynamic  efficiency  defined  as  the  final  kinetic  energy  of  the  foil  divided  by  the  total  energy  absorbed  by 
the  target,  U  \ s  the  average  ablation  velocity,  and  ^is  the  average  target  velocity. 

Both  solutions  A  and  B  are  initialized  at  an  absorbed  laser  intensity  of  8  x  1011  w/cm2, 
corresponding  to  a  time  of  3.2  nsec  before  the  peak  intensity.  We  model  the  time  evolution  of  the 
laser  pulse  with  a  Gaussian  having  the  appropriate  width.  These  two  initial  quasi-static  solutions  were 
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chosen  as  spanning  the  kind  of  different  steady  state  results  admissible  in  our  profile  generating  equa¬ 
tions.8  This  ambiguity  arises  irom  such  things  as  the  choice  of  the  adiabat  of  the  foil  and  the  details  of 
the  increase  of  the  laser  intensity.  The  initial  density  profiles  for  these  two  solutions  are  shown  in  Fig. 
B2b.l. 

The  first  four  rows  of  Table  B2b.I  display  the  results  obtained  using  the  ID  code.  Comparing  the 
first  and  third  rows  with  each  other  and  with  the  experimental  results  of  the  last  row,  we  see  that  the 
differences  between  solution  A  and  B  are  minimal,  and,  also,  that  the  agreement  with  the  experimental 
mass  loss  and  ablation  velocity  is  good.  However,  both  the  target  velocity  and  hydrodynamic  efficiency 
are  calculated  to  be  too  large.  This  discrepancy  arises  since  we  are  performing  a  ID  calculation.  In  the 
calculation,  all  the  momentum  of  the  ablated  plasma  is  imparted  to  the  foil.  In  the  experiment,  how¬ 
ever,  the  ablated  material  is  seen  to  emerge  in  a  cone  with  an  opening  angle  of  roughly  20°.  Hence,  we 
are  overestimating  the  impulse  to  the  target  that  the  ablated  matter  transmits  and  our  calculated  target 
velocity  and,  hence,  efficiencies  are  too  high.  Results  of  our  simplified  2D  calculations  will  be  dis¬ 
cussed  shortly. 

For  the  ID  calculations,  runs  where  the  classical  plasma  thermal  conductivity  coefficient  was  arbi¬ 
trarily  reduced  everywhere  by  a  factor  of  ten  were  performed.  These  results  are  displayed  in  rows  two 
and  four.  Somewhat  surprisingly,  for  a  factor  of  ten  change  in  the  thermal  conductivity,  only  relatively 
minor  changes  occur  in  the  results.  As  expected,  the  ablation  region  becomes  hotter  and  the  ablation 
velocity  increases  for  the  lower  thermal  conductivity  case.  Also,  as  expected  for  this  case,  less  mass  is 
ablated,  the  ablation  is  less  efficient,  and  the  target  achieves  a  lower  velocity. 

Part  of  the  reason  for  this  small  sensitivity  to  the  thermal  conduction  coefficient  can  be  deduced 
from  Fig.  B2b.2.  Here  for  solution  A,  the  density  profiles  at  the  time  of  peak  intensity  are  displayed  for 
the  two  different  values  of  thermal  conductivity.  For  the  lower  value,  the  density  gradients  have 
dramatically  steepened  so  that  the  distance  from  the  critical  density  region  to  the  peak  density  is  much 
shorter.  Hence,  the  mass  contained  between  peak  and  critical  density  is  greatly  reduced  and  this,  to  a 

large  extent,  compensates  for  the  reduced  conductivity. 
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As  a  test  of  the  consistency  of  our  two  codes,  we  performed  a  calculation  with  the  2D  code  in  a 
ID  mode.  All  the  quantities  were  initialized  with  no  ^-variation  (the  w-direction  being  transverse  to  the 
laser  beam)  and  this  uniformity  persisted  throughout  the  run.  The  results  of  this  run  are  displayed  in 
the  fifth  row  of  Table  B2b.I.  The  agreement  between  the  first  and  fifth  rows  is  excellent. 

Next,  we  simulated  2D  effects  by  adding  a  v-variation  to  the  laser  intensity.  The  intensity  was 
uniform  in  a  region  240  n  wide  (approximately  the  experimental  spot  diameter)  and  then  fell  off 
exponentially  on  either  side.  It  was  much  more  difficult  to  extract  the  numbers  needed  for  a  compari¬ 
son  to  experiment  in  this  case,  and  only  the  entries  for  the  averaged  quantities,  the  two  velocities,  have 
been  filled  in  for  the  sixth  row  of  Table  B2b.I.  It  was  too  difficult  to  accurately  locate  all  of  the  foil 
mass  and  so  the  summed  quantities  were  not  calculated.  However,  the  two  velocities  are  sufficient  for 
our  purposes.  Note  that  there  is  a  significant  drop  in  the  target  velocity,  as  expected,  when  2D  effects 
are  included.  The  fact  that  there  was  such  a  small  change  in  the  ablation  velocity,  again  as  expected, 
gives  us  confidence  that  we  were  able  to  extract  these  average  velocities  correctly  from  our  2D  calcula¬ 
tion. 


All  of  these  results  were  extracted  from  the  computed  results  at  a  time  corresponding  to  about  5.5 
nsec  after  the  peak  intensity  was  reached.  At  that  time,  the  laser  intensity  was  down  about  a  factor  of 
fifty  and  only  minor  changes  should  occur  at  later  times.  That  this  was  indeed  the  case  was  verified  in  a 
ID  run  that  was  extended  to  18  nsec  beyond  the  intensity  peak  (corresponding  to  a  flux  reduction  of 
over  35  orders  of  magnitude). 

We  have  found  that  the  internal  consistency  between  a  ID  and  a  2D  computer  model  to  calculate 
the  ablative  acceleration  of  a  thin  foil  by  a  relatively  long  pulse,  low  intensity  laser  beam  is  excellent. 
A  purely  ID  calculation  can  adequately  describe  the  mass  burn-off  and  blow-off  velocity,  but  2D  effects 
or  a  variable  metric  ID  geometry  must  be  included  to  calculate  the  target  velocity  to  better  than  30% 


accuracy. 
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As  long  as  reasonable  initial  profiles  are  used  and  the  thermal  conduction  is  within  an  order  of 
magnitude  of  the  classical  value,  it  appears  that  conservation  of  mass,  momentum  and  energy  are 
sufficient  to  give  reasonable  agreement  with  experiment.  However,  detailed  understanding  of  the 
microscopic  processes  at  work  requires  more  accurate  measurements  and  modeling  as  well  as  more 
information  about  the  density  and  temperature  of  the  ablation  region. 

REFERENCES 

1.  NRL  Laser-Plasma  Interaction  Group,  NRL  Memo  Report  No.  3890  (1978),  ed.  B.H.  Ripin;  R. 

Decoste  et  aL,  Phys.  Rev.  Lett.  42,  1673  (1979);  B.H.  Ripin  et  al.,  Phys.  Rev.  Lett.  43,  350 

(1979). 

2.  O.N.  Krokhin  et.  al.,  Sov.  Phys.  JETP  42,  107  (1976);  N.G.  Basov  et.  al.,  Sov.  Phys.  JETP  Lett. 

26,  433  (1977);  J.P.  Anthes,  M.A.  Gusinow,  and  M.  Keith  Matzen,  Phys.  Rev.  Lett.  41,  1300 

j  (1978). 

I 

j 

3.  J.M.  McMahon,  NRL  Memo  Report  No.  3591  (1977),  ed.  S.E.  Bodner. 

4.  R.E.  Kidder,  Nucl.  Fusion  8,  3  (1968);  K.A.  Brueckner,  P.M.  Campbell,  and  R.A.  Grandy,  Nucl. 
Fusion  15,  471  (1975);  S.J.  Gitomer,  R.L.  Morse,  and  B.S.  Newberger,  Phys.  Fluids  20,  234 
(1977);  R.C.  McCrory  and  R.L.  Morse,  Phys.  Rev.  Lett.  38,  544  (1977);  F.S.  Felber,  Phys.  Rev. 
Letter,  39,  84  (1977);  J.D.  Perez  and  G.L.  Payne,  Phys.  Fluids  22,  361  (1979). 

5.  J.P.  Boris,  NRL  Memo  Report  No.  4022  (1979). 

6.  J.P.  Boris,  and  D.L.  Book  in  "Methods  of  Computational  Physics"  16,  85  (1976);  J.P.  Boris  in 
"Comments  on  Plasma  Physics  and  Controlled  Fusion"  3,  1  (1977);  Inertial  Confinement  Fusion 

i  at  NRL,  S.  Bodner,  et  al.,  IAEA-CN-37-B-3  Proceedings  of  the  Innsbruck  IAEA  Conference, 

1978. 

7.  S.I.  Braginskii  in  "Review  of  Plasma  Physics"  1,  205,  ed.  M.A.  Leontovich. 

8.  J.H.  Orens,  NRL  Memo  Report  4167  (1980). 


40 


NRL  MEMORANDUM  REPORT  4369 


Table  B2b.I 

A  Comparison  of  Calculations  to  Experiment  for  a  15#* 
CH  Foil  Irradiated  with  a  1013  w/cm2  Peak  Intensity 
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Fig.  B2b.I— Two  very  different  initial  density  profiles  corresponding  to  an  absorbed  laser 
intensity  of  8  x  1011  w/cm2  and  a  15m  CH  target.  This  is  incident  from  the  right. 
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Fig.  B2b  2-Two  density  profiles  at  the  time  of  the  peak  intensity  that  both  evolved  from  initial  solution  A.  Only  the  thermal 
conduction  coefficients  differ  for  the  two  cases  as  labelled.  Note  the  large  difference  in  the  location  of  the  critical  density  relative 
to  the  density  peak.  This  figure  is  based  on  results  from  our  ID  calculations.  Again,  the  laser  is  incident  from  the  right. 
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C  DEVELOPMENT  OF  COMPUTATIONAL  TECHNIQUES 

Realistic  numerical  calculations  of  the  Rayleigh-Taylor  instability  at  interior  material  interfaces 
require  accurate  advancement  of  the  material  interfaces.  This  requirement  favors  the  use  of  a  Lagran- 
gian  formulation  because  Eulerian  techniques  must  either  employ  special  procedures  to  locate  an  inter¬ 
face  or  use  excessively  fine  gridding  to  resolve  it.  Unfortunately  these  Eulerian  procedures  are  not  only 
cumbersome  and  time  consuming,  but  result  in  a  numerical  diffusion  which  can  grow  large  enough  to 
mask  the  physical  motion  of  the  interface.  The  price  that  is  paid  for  ease  in  locating  an  interface  in  a 
Lagrangian  treatment  is  a  distorted  or  skewed  grid  in  all  but  the  simplest  of  flows. 

A  badly  distorted  grid  can  always  be  regularized  by  an  Eulerian  rezone  phase,  but  only  at  the 
expense  of  reintroducing  numerical  diffusion.  The  accuracy  of  these  calculation  is  then  limited  to  the 
accuracy  of  the  rezone  algorithms,  which  are  often  less  accurate  than  the  corresponding  strictly  Eulerian 
treatments.  Since  the  purpose  of  rezoning  is  to  maintain  a  roughly  central  position  for  a  vertex  relative 
to  its  neighbors,  an  alternative  solution  would  be  to  allow  the  grid  to  reconnect  as  vertices  migrate  in 
the  flow  to  reflect  their  new  relative  positions.  Unfortunately,  a  single  line  reconnection  on  a  quadrila¬ 
teral  grid  replaces  two  quadrilateral  cells  with  one  triangular  and  one  pentagonal  cell,  an  therefore 
requires  special  differencing  algorithm  in  the  neighborhood  of  the  reconnection.  In  general  this  means 
that  grid  reconnections  for  a  quadrilateral  mesh  must  be  global  and  can  result  in  a  quadrilateral  mesh 
only  for  periodic  grid  anomalies. 

in  contrast,  if  a  triangular  grid  is  used  instead  of  a  quadrilateral  mesh,  any  reconnected  grid  line 
must  result  in  two  triangles.  Grid  anomalies  can  be  resolved  locally  without  recourse  to  global  rezoning 
and  cannot  result  in  changing  the  liaracter  of  grid  zones.  Vertices  can  also  be  added  or  deleted  where 
required  to  maintain  resolution  with  only  local  grid  changes.  In  other  words,  it  is  possible  to  dynami¬ 
cally  alter  a  triangular  grid  at  a  given  location  with  the  effects  of  the  alteration  being  felt  only  by  the 
vertices  near  that  location.  This  fact  opens  up  the  possibility  of  conserving  physical  variables  locally 
even  during  grid  restructuring.  In  the  hydrodynamics  code  SPLISH  this  goal  is  attained  through  the  use 
of  conservation  laws  expressed  as  integrals  over  regions  about  each  vertex  and  over  triangles  or  pairs  of 
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triangles.  Grid  reconnections,  vertex  additions  or  vertex  deletions  result  in  the  unambigious  definition 
of  new  cell  quantities  through  the  use  of  these  integral  invariants.  Regridding  algorithms  can  therefore 
be  formulated  which  are  completely  automated  and  local  and  which  preserve  conserved  physical  vari¬ 
ables. 

The  incompressible  fluid  code  SPLISH  incorporates  these  algorithms  into  a  general  connectivity 
hydrodynamics  code  using  a  triangular  grid.  Calculations  of  the  Rayleigh-Taylor  instability  with  SPLISH 
had  previously  been  extended  well  into  the  nonlinear  regime.  Figure  C.l  illustrates  one  earlier  calcula¬ 
tion  which  modeled  the  effect  of  the  unstable  pusher-ablator  interface  on  the  stable  pusher-fuel  inter¬ 
face.  The  lowest  region  in  the  figure  represents  the  ablator  which  is  separated  from  the  pusher  by  an 
interface  that  is  given  an  initial  sinusoidal  perturbation.  The  fuel  is  represented  by  the  ungridded 
region  above  the  pusher  and  is  separated  from  it  by  a  free  surface.  The  shape  of  the  interface  at  four 
times  in  the  calculation  is  shown  in  both  the  plots  of  the  triangular  grid  and  the  pathlines.  Pathline 
plots  record  the  current  and  five  previous  positions  of  each  vertex.  The  development  of  Kelvin- 
Helmholtz  vortices  on  the  interface  due  to  the  progressively  strengthening  shear  is  evident  in  the  path¬ 
lines. 

The  addition  of  gridding  algorithms  to  permit  attachment  of  the  heavier  fluid  at  the  bottom  of  the 
computational  region  extended  these  runs  considerably.  The  accelerating  columns  of  dense  pusher 
eventually  replaced  the  ablator  along  most  of  the  ablation  front,  leaving  a  thin  column  of  ablator  being 
pumped  toward  the  fuel  by  the  vortex  pair.  A  complex  flow  of  mixed  pusher  and  ablator  continued  to 
alter  the  vorticity  "baroclinically"  at  each  of  the  vertices  and  a  jet  of  pusher  was  shown  to  impinge  into 
the  fuel  region.  At  this  point  in  the  calculation  higher  resolution  was  desirable  as  well  as  an  improved 
treatment  of  reattachment  at  the  free  surface.  Therefore  the  fuel  region  itself  was  gridded  and  runs 
were  begun  with  better  resolution. 

Gridding  beyond  the  free  surface  and  within  cavities  simplifies  the  grid  restructuring  problems 
there  but  leaves  unresolved  the  question  of  improved  boundary  conditions  at  the  free  surface.  At  a 
free  surface  pressures  are  assumed  to  be  given.  Although  this  defines  pressure  gradients  as  well,  the 


LABORATORY  FOR  COMPUTATIONAL  PHYSICS  AND  PLASMA  PHYSICS  DIVISION 

boundary  condition  is  not  yet  sufficiently  specified.  Pressures  are  used  in  the  interior  of  the  fluid  to 
maintain  a  flow  consistent  with  the  chosen  equation  of  state  (V  •  V  —  0  in  SPL1SH).  Therefore,  at  the 
surface  some  other  means  of  satisfying  this  condition  must  be  found  for  the  partial  surface  vertex  cells. 

Currently  the  code  adjusts  the  free  surface  positions  to  keep  the  surface  layer  divergence-free. 
Therefore  the  free  surface  boundary  in  now  treated  in  a  manner  analogous  to  that  used  at  rigid 
boundaries— by  vertex  movement.  Vertices  at  a  rigid  boundary  are  moved  to  locations  enforced  by  the 
contour  of  the  boundary.  At  the  free  surfaces,  movement  of  the  vertices  provides  for  divergence-free 
flow  about  the  vertex.  This  boundary  condition  has  now  been  implemented  and  has  eliminated  a  previ¬ 
ous  mode  of  failure  due  to  the  presence  of  tangential  velocities  at  the  free  surface.  Although  this  for¬ 
mulation  is  more  consistent  and  improves  the  free  surface  behavior,  the  free  surface  motion  is  still 
problematic  in  that  a  first  order  approximation  is  used.  Higher  order  solutions  are  currently  being 
sought. 

Improved  resolution  for  the  calculations  was  also  needed,  but  could  be  achieved  only  by  recoding. 
Although  the  number  of  vertices  could  be  increased  easily  by  a  factor  of  4,  this  is  only  a  doubling  of 
resolution  in  either  direction.  The  construction  of  a  larger  code  evolved  into  three  major  work  areas. 
First,  large  outdated  portions  of  old  code  were  cleaned  out.  This  includes  the  old  method  of  marker 
particles  which  had  been  superceded  by  the  current  pathlines  as  well  as  the  elimination  of  large 
amounts  of  common  area  which  was  being  used  inefficiently.  The  decrease  in  core  needed  for  the  pro¬ 
gram  increased  the  total  number  of  vertices  that  can  be  used  in  a  practical  calculation  to  roughly  2000, 
almost  50%  more  than  possible  before  the  clean-up. 

As  Figure  C.2  illustrates,  the  finer  resolution  runs  permitted  much  improved  detail  for  several 
features  of  the  late-time  flow.  A  thin  ribbon  of  pusher  remains  between  the  ablator  and  fuel  at  this 
time.  The  jet  formed  above  the  main  column  of  pusher  has  collapsed,  entraining  two  bubbles  of  fuel. 
Similar  "bubbles"  of  pusher  within  the  ablator  and  ablator  within  the  pusher  are  found  within  the  vortex 
pair.  The  pathlines  are  of  course  much  shorter  in  this  simulation  because  of  the  decreased  timestep 
necessary  at  this  point  in  the  calculation  due  to  the  Courant  condition. 
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Figures  C.3  and  C.4  illustrate  another  calculation  in  which  an  initial  mode  2  perturbation  to  the 
interface  contained  a  slight  mode  3  admixture.  The  gross  behavior  of  the  interface  is  similar  to  the  ear¬ 
lier  calculations  except  that  the  two  '’spikes”  no  longer  progress  at  the  same  rate.  Furthermore,  it  is  evi¬ 
dent  in  the  last  frame  that  the  higher  mode  structure  has  accelerated  the  growth  of  secondary  vortices 
near  the  thinned  pusher. 

Second,  work  has  been  undertaken  on  improving  the  elliptic  solver.  Particularly  at  large  array 
sizes,  the  elliptic  solver  overwhelmingly  dominates  the  cost  of  a  simulation,  and  a  more  rapid  conver¬ 
gence  could  make  very  large  array  sizes  practical.  As  a  first  attempt  at  such  an  improvement,  a  totallv 
scalar  version  was  implemented  which  used  SOR  and  which  improved  the  rate  of  convergence  substan¬ 
tially.  However,  after  several  tests  this  approach  was  abandoned,  since  the  gains  in  rate  of  convergence 
were  offset  by  the  greatly  increased  computation  times.  A  vectorized  version  executes  much  more 
rapidly  on  the  TI-ASC,  so  that  the  vectorized,  poorly  converging  algorithm  is  far  less  costly  than  the 
scalar,  rapidly  converging  one.  Efforts  are  now  underway  to  find  other  means  of  solving  the  elliptic 
equations  which  are  more  adapted  to  the  particular  problems  presented  by  a  general  triangular  mesh. 

Although  our  results  indicate  that  the  code  is  nearly  ready  for  a  systematic  study  of  the  influence 
of  the  Rayleigh-Taylor  instability  on  the  ablation  dynamics,  some  code  modifications  still  need  to  be 
performed.  The  analysis  of  the  pressure  fields  generated  in  current  runs  indicates  fluctuation  due  to 
variations  in  the  velocity  fields.  These  velocity  fluctuations  have  in  turn  been  traced  to  the  rotation 
operator  used  in  SPL1SH  to  conserve  vorticity.  Exact  conservation  of  vorticity  in  a  Lagrangian  code 
was  attempted  for  the  first  time  in  SPLISH.  This  exact  conservation  narrows  the  choice  of  rotation 
operators  to  a  single  possible  operator.  However,  this  exact  operator  can  induce  the  pressure  fluctua¬ 
tions  described  above  in  regions  of  large  shear.  Vertices  in  the  shear  layer  experience  an  added  diver¬ 
gence  from  the  rotation  operator.  Since  this  added  divergence  is  numerical  and  not  physical,  SPLISH  is 
currently  being  altered  to  iteratively  reduce  these  added  divergences.  Such  a  solution  retains  the  exact 
conservation  of  vorticity  yet  eliminates  any  spurious  numerical  divergence  in  regions  of  shear. 
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While  this  modification  is  being  implemented,  another  alteration  is  also  being  coded.  Velocities 
are  placed  at  the  triangle  centroids  in  SPLISH.  During  reconnection,  a  quadrilateral  velocity  is  con¬ 
served,  as  well  as  vorticity  and  divergence  about  all  vertices.  However,  despite  this  conservation  each 
of  the  four  affected  vertices  experience  a  net  acceleration  due  to  the  weighting  procedure  used  in 
assigning  vertex  velocities.  A  new  weighting  procedure  is  being  implemented  concurrently  with  the 
new  rotation  operator  to  eliminate  these  accelerations. 

Both  of  these  alterations  were  found  to  be  necessary  primarily  through  the  use  of  new  diagnostic 
routines  for  the  velocity  fields.  These  routines  were  formulated  to  aid  in  debugging  the  added  gridding 
routines  mentioned  previously.  A  side  effect  of  their  use  should  be  improved  convergence  for  our 
elliptic  solver  due  to  the  elimination  of  divergences  introduced  through  both  the  rotator  and  velocity 
weighting.  The  addition  of  these  divergences  to  the  physical  flow-induced  divergences  implies  spurious 
pressure  changes  which  could  be  eliminated  previously  only  by  more  iterations  in  the  solver. 

Efforts  are  also  underway  to  construct  a  fully  single  precision  version  of  the  code  to  both  decrease 
the  CPU  time  and  storage  currently  used  by  SPLISH.  Although  the  bulk  of  the  gridding  routines  are 
currently  single  precision,  the  main  hydrodynamic  routines  have  been  double  precision.  This  had 
mistakenly  been  thought  necessary  in  the  past  to  ensure  convergence  of  the  elliptic  solver.  Preliminary 
tests  indicate  that  the  new  rotation  operator  eliminates  the  need  for  double  precision  in  these  routines. 

And  third,  an  effort  has  been  initiated  to  standardize  the  routines  used  in  this  project.  A  change 
of  array  sizes  or  particular  problem  generally  has  meant  a  change  to  most  subroutines  in  the  code.  This 
led  to  multiple  versions  of  source  code.  A  perhaps  inevitable  result  of  multiple  versions  has  been  their 
gradual  evolution  and  mutation  into  different  codes  instead  of  different  versions  of  the  same  code,  and 
an  accompanying  rise  in  errors  due  to  this  multiplicity.  Standardization  will  eliminate  this  source  of 
error 

Several  presentations  have  been  given  on  the  results  of  the  Lagrangian  calculations.  These 
include  the  following  oral  papers: 
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"Simulations  of  the  Nonlinear  Development  of  the  Rayieigh-Taylor  Instability  for  Iner¬ 
tial  Confinement  Fusion"  at  the  Boston  APS  meeting,  1979. 

"Occurrence  of  the  Rayieigh-Taylor  Instability  in  Inertial  Confinement  Fusion"  at  the 
Fluid  Dynamics  sectional  APS  meeting  at  Notre  Dame,  1979. 
and 

"Some  Considerations  of  Rayieigh-Taylor  Instability  in  Inertial  Confinement  Fusion"  at 
the  1980  CLEOS/1CF  meeting. 

In  addition,  the  calculations  discussed  here  are  featured  in  a  chapter  on  the  triangular  grid  technique  in 
a  forthcoming  book.  Recent  Developments  in  Computational  Techniques  for  Applied  Hydrodynamics , 
Springer-Verlag. 
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Fig.  C.l  —  Flow  patterns,  interface  location  and  grid  structure  of  the  internal  interface  Rayleigh-Taylor  instability  as  computed  by 
the  SPLISH  computer  code  The  upper  left  quadrant  of  the  figure  shows  the  initial  grid  of  triangles  with  a  small  sinusoidal  per¬ 
turbation  of  the  interface  at  t  —  0  The  initial  velocities  (in  the  accelerating  frame  of  reference)  were  all  taken  to  be  zero.  The 
lower  left  quadrant  shows  the  flow  field  as  the  instability  enters  its  nonlinear  phase.  The  vertex  locations  are  plotted  at  successive 
timesteps  to  give  streaks  whose  length  is  proportional  to  the  local  fluid  velocity.  Shear  flow  is  established  as  the  interface 
steepens,  giving  rise  to  vortices  centered  on  the  interface  in  a  flow  which  has  been  likened  to  the  onset  of  the  Kelvin-Helmholtz 
instability  In  the  two  right-hand  quadrants,  two  fluids  are  beginning  to  mix  in  the  Rayleigh-Taylor  vortices  and  the  rate  of  fall  of 
the  heavy  material  is  no  longer  increasing. 
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g.  C.2  —  The  Rayleigh-Taylor  instability  at  a  later  time  using  finer  resolution.  Evidence  of  smaller  scale  structure  includes 
secondary  vortices,  collapsing  jets,  waves  on  the  ribbon  of  pusher,  and  entrained  "bubbles"  of  fuel,  ablator,  and  pusher. 
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dominant  mode  2  perturbation.  An  admixture  of  mode  3  is  evident  in  the 
asymmetry  in  both  the  .Y  and  Y  directions. 
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0.0  l  .  0 

Fig.  C.4  —  Nonlinear  development  of  the  mixed  mode  2  and  mode  3  calculation  of 
Fig.  C.3.  The  mode  3  admixture  has  led  to  different  levels  of  entrainment  for  each 
spike,  but  the  vertical  descent  is  roughly  similar. 
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Dl.  Discussion  of  the  Quasi-Static  Equilibrium  Parameter  Space 

In  our  earlier  work1;  we  developed  analytic  approximations  and  numerical  solutions  for  the  struc¬ 
ture  of  quasi-static  laser  driven  ablation  layers.  At  that  time  we  discussed  the  dependence  of  the  struc¬ 
ture  on  two  independent  dimensionless  parameters  M0  and  iV0.  Af0  is  the  isothermal  mach  number  at 
the  density  peak,  which  for  our  situation  is  usually  a  small  number,  and  N0  is  the  Peclet  number  at  the 
density  peak  (thermal  diffusion  time  scale/free  streaming  time  scale),  which  in  our  regime  is  usually  a 
large  number.  Since  thermal  diffusion  and  free  streaming  are  dominant  mechanisms  in  our  ablation 
model  it  is  reasonable  they  should  parameterize  the  solution.  Although  A/0  and  JV0  are  the  most  con¬ 
venient  parameters  for  the  analysis  it  will  be  shown  subsequently  that  the  acceleration  and  entropy  of 
the  shell  can  equally  well  parameterize  the  solution. 

Within  certain  limitations,  the  choice  of  values  for  Af0  and  /V0  is  arbitrary.  In  essence,  there  is 
not  a  unique  equilibrium  for  the  prescribed  physical  constraints,  laser  intensity,  and  shell  mass,  but 
various  accessible  equilibria.  Exactly  which  equilibria  would  occur  for  a  given  physical  situation  might 
be  determined  by  the  necessarily  time-dependent  evolution  leading  up  to  that  equilibrium.  Two  such 
effects  might  be  the  entropy  and  acceleration  of  the  shell  imparted  by  an  initial  shock.  Recently  we 
analyzed  the  A/0  —  N0  parameter  space  for  curves  of  constant  shell  entropy  and  acceleration.  The 
results  of  such  an  analysis  for  15  ^  C-H  and  an  intensity  of  1013  w/cm2  are  given  in  Fig.  Dl.l.  The 
curves  of  constant  entropy  are  nearly  straight  lines  with  slopes  approximately  inversely  proportional  to 
the  magnitude  of  the  entropy.  Initially,  visual  inspection  seemed  to  imply  that  all  the  straight  lines 
emanated  from  a  point  on  the  negative  A/0  axis,  but  more  detailed  analysis  seems  to  show  that  the 
curves  are  not  straight  near  A/0  -  0  and  begin  to  bend  toward  the  origin.  Further  analysis  will  be 
necessary  before  the  exact  intersection  point  can  be  determined.  The  curves  of  constant  acceleration 
are  somewhat  like  hyperbolas,  with  the  curves  for  increasing  acceleration  being  further  from  the  origin. 
Note  that  only  a  small  variation  in  the  acceleration  causes  large  variations  in  M0  and  iV0.  In  our  earlier 
work  we  identified  two  regions  for  the  pressure  profile.  Below  the  curve  A0A/q  -  .4\/3  (shown  as  the 
lower  dashed  curve  on  Fig.  Dl.l)  the  pressure  profile  cannot  have  a  turning  point  and  the  pressure 
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increases  monotonically  out  to  the  critical  surface.  This  would  correspond  to  a  case  where  the  density 
and  pressure  gradients  are  always  opposite  in  sign  from  the  density  maximum  until  the  critical  surface. 
Such  a  case  would  be  expected  to  be  Rayleigh-Taylor  unstable.  Above  this  curve  the  pressure  profile 
can  turn  over  and  yield  a  region  where  the  density  and  pressure  gradients  are  again  of  like  sign.  Recent 
investigation  has  shown  that  there  is  a  transition  region  between  the  case  where  the  pressure  only 
increases  and  where  it  reaches  a  maximum  and  then  decreases  out  to  the  critical  surface.  Previously  we 
have  shown  that 
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Therefore  in  the  former  case  — never  gets  as  large  as  — 77  while  in  the  latter  case  once  it  does 
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become  larger  it  never  becomes  smaller  again  before  reaching  the  critical  surface.  Recent  numerical 

analysis  has  observed  the  transition  between  the  two  cases  where  does  again  become  smaller  than 

«{ 

—77  prior  to  the  critical  surface  and  the  pressure  profile  again  increases.  This  region  is  enclosed 
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between  the  two  dashed  curves  on  Fig.  Dl.l.  As  the  equilibrium  moves  from  the  lower  dashed  curve 
to  the  upper  one,  the  second  turning  point  moves  closer  to  the  critical  surface.  Therefore  by  moving 
from  the  lower  to  the  upper  dashed  curve  the  overall  region  where  the  density  and  pressure  gradients 
are  opposite  in  sign  becomes  smaller.  Above  the  upper  dashed  curve  there  is  only  a  very  small  region 
of  opposed  gradients  and  this  is  the  profile  considered  in  most  previous  analyses.  It  is  interesting  to 
note  that  the  shape  of  these  dashed  curves  very  much  resembles  the  curves  for  constant  acceleration. 
Therefore  the  results  seem  to  imply  that  below  a  certain  acceleration  the  pressure  never  turns  over 
while  above  a  certain  acceleration  there  is  only  one  turning  point  for  the  pressure.  Of  course,  further 
study  is  required  before  making  a  firm  statement.  Outside  the  largest  constant  acceleration  curve  on 
Fig.  Dl.l,  the  flow  out  of  the  critical  surface  is  no  longer  supersonic.  This  region  is  excluded  from  the 
analysis  since  subsonic  flow  in  the  underdense  gas  would  require  the  temperature  to  increase  away  from 
the  heat  source  at  xc  which  would  be  unphysical.2  Again  the  onset  of  this  region  seems  to  be  strongly 
dependent  on  the  acceleration. 
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Since  the  determination  of  a  specific  equilibrium  depends  on  the  entropy  and  acceleration  of  the 
foil  which  are  not,  in  general,  specifiable  quantities  we  also  investigated  the  possible  existence  of  evolu¬ 
tionary  trends  for  the  various  types  of  "quasi-static”  pressure  profiles.  Using  idealized  fully  time- 
dependent  calculations  at  constant  laser  intensity  we  have  noted  that  profiles  in  which  the  pressure 
turns  over  seem  to  retain  their  basic  form,  while  profiles  in  which  the  pressure  continually  increases 
seem  to  evolve  into  profiles  in  which  the  pressure  turns  over.  We  conjecture  that  the  latter  solution 
evolves  such  that  the  entropy  of  the  foil  remains  constant  while  the  acceleration  of  the  foil  increases 
into  a  region  where  the  pressure  would  then  turn  over.  Figure  D1.2a  shows  the  temporal  evolution  of 
a  "quasi-static”  equilibrium  with  a  continually  increasing  pressure  profile,  while  Fig.  D1.2b  shows  the 
relative  invariance  of  a  "quasi-static”  equilibrium  with  a  pressure  profile  that  turns  over.  These  evolu¬ 
tionary  trends  are  very  interesting  since  they  seem  to  imply  that  the  only  physically  realizable  equilibria 
might  be  those  in  which  the  pressure  profile  has  a  turning  point  and  hence  is  most  stable  to  the 
Rayleigh-Taylor  modes. 
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Fig.  DM  —  Contours  for  constant  entropy  (nearly  straight  lines)  and  constant  acceleration  (curved  lines)  in  A#o~.V0  space  for 
the  case  of  a  15^  plastic  foil  irradiated  with  tO13  w/cm:  laser  intensity  (80%  absorption).  Entropy  contours  are  labeled  in  units  of 
10,J  and  acceleration  contours  are  labeled  in  units  of  lO^cm/sec2.  Dashed  lines  approximately  delineate  the  regions  of  stable  and 
unstable  cases. 
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D2.  The  Vorticity  Generation  Model 

One  of  the  most  important  areas  in  inertial  confinement  laser  fusion  research  concerns  the 
Rayleigh-Taylor  stability  of  the  ablation  layer  for  laser  driven  foils  and  spherical  implosions.  Since  the 
density  and  pressure  gradients  are  opposite  in  sign  for  that  region,  simple  analysis  predicts  that  pertur¬ 
bations  should  grow  exponentially  there.  Therefore  it  is  of  great  interest  to  determine  under  what  cir¬ 
cumstances  the  expected  instability  can  be  stabilized.  Our  earlier  research  into  the  quasi-static  (time 
independent)  ablation  layer  profiles  has  identified  two  distinct  types  of  pressure  profiles.  In  the  first 
case  the  pressure  has  the  intuitively  expected  turning  point  and  yields  a  very  narrow  region  where  the 
density  and  pressure  gradients  are  opposite  in  sign.  In  the  second  case  the  pressure  does  not  turn  over 
and  continues  to  increase  out  to  the  critical  surface.  Here  the  region  where  the  density  and  pressure 
gradients  are  opposite  in  sign  is  quite  broad.  There  also  is  a  transition  region  where  the  pressure  does 
have  a  local  maximum  near  the  density  peak,  but  then  at  a  point  between  there  and  the  critical  surface 
the  pressure  increases  again.  Recently  we  have  utilized  a  vorticity  generation  formulation  to  study  the 
stability  of  these  cases  and  preliminary  results  show  that  the  first  case  is  stable  while  the  second  is 
unstable.  Complex  interactions  among  the  various  external  parameters  determine  whether  or  not  the 
pressure  profile  turns  over,  but  the  predominant  indicator  is  the  maximum  magnitude  of  the  slope  in 
the  ablation  region  of  the  dimensionless  density  profile. 

In  the  usual  analysis  the  vorticity  equation  is  formed  by  taking  the  curl  of  the  momentum  equa¬ 
tion  to  obtain 

+  V  •  (V£)  -  (  VV  -  (1) 

Ot  pL 

where  f  =  VxV  is  the  vorticity,  p  is  the  density,  P  is  the  pressure,  and  the  source  term  on  the  right 
hand  side  is  the  baroclinic  vorticity  generation  term.  In  two  dimensions,  here  treated  as  cartesian,  f 
reduces  to  a  single  scalar  component  which  will  be  taken  along  the  z  axis.  Thus  the  third  term  on  the 
left  hand  side  of  Eq.  (1)  vanishes  and  the  equation  reduces  to: 
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|f  +  V  •  (Vf )  -  ^P*V-  (2) 

of 

where  (Ax,  y,  /)  =  (Vx  V)z  -  and  u  and  v  are  respectively  velocities  in  the  x  and  y  direc¬ 

tions.  Linearizing  Eq.  (2)  with  respect  to  an  equilibruim  flow  only  in  the  x  direction  and  independent 


of  y  and  z  yields  for  the  first  order  perturbation  (tilde)  quantities 
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Since  for  the  equilibrium 


(3) 


Eq.  (3)  may  be  written 
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After  some  manipulation  and  using  the  first  order  y  momentum  equation 


Eq.  (7)  becomes 
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It  is  interesting  to  note  that  the  pressure  perturbation  does  not  appear  in  Eq.  (7). 


(7) 


In  studies  of  Rayleigh-Taylor  stability  it  is  usually  sufficient  to  consider  incompressible  perturba¬ 
tions.  Furth,  Killeen  and  Rosenbluth3  and  Coppi,  Green  and  Johnson4  used  this  approach  for  plasma 
slab  stability  and  Kidder  justified  it  for  isentropic  implosions.5  Thus  the  first  order  velocity  can  be 
derived  from  the  curl  of  a  vector  potential. 

S  =  £  andis-it.  ,8) 

dy  dx 


Letting 
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yields  a  potential  equation  of  the  form 
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Differentiating  the  first  order  continuity  equation 


with  respect  to  y  yields 
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dr 

where 

In  this  representation  Eq.  (7)  becomes 
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Considering  perturbations  of  the  form 


(12) 


(13) 


/»/  (x,  r )  exp  Uky).  (14) 

Eqs.  (10),  (12),  and  (13)  become 


(15a) 

(15b) 

1 2  9«  r 

-*<>  »?*• 

(15c) 

Equations  (15a),  (15b),  and  (15c)  become  the  basis  for  the  vorticity  generation  model.  Nowhere  in  the 
analysis  have  the  zero  order  quantities  been  assumed  to  be  independent  of  time  so  these  equations  can 
just  as  well  be  applied  to  study  the  stability  of  time  dependent  situations.  If  the  condition  of  an 
incompressible  perturbation  is  relaxed  then  Eq.  (15a)  would  take  the  form  of  a  wave  equation  allowing 
the  presence  of  sound  waves  to  govern  the  onset  of  motion  at  one  point  after  vorticity  is  generated  at 


another. 
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D2a.  Vorticity  Generation  Model  Tests  and  Application  to  Laser  Driven  Ablation  Profiles 


Prior  to  applying  our  vorticity  generation  model  to  the  actual  laser  ablation  profiles  we  tested  the 
method  in  cases  of  an  accelerating  step  density  discontinuity  where  the  answers  are  known  analytically. 
Figure  D2a.l  shows  the  density  profile  for  the  cases  chosen.  Equilibria  both  with  and  without  mass 
flow  were  considered.  In  the  cases  without  mass  flow  the  growth  rate  of  the  Rayleigh-Taylor  modes  is 
given  by  the  standard  formula 

y  ■*  \TkgA  (16a) 

where 


Previous  authors6  7 


A  =  P-C 

Pa  +  Pb 

have  considered  the  cases  with  constant  mass  flow  and 


(16b) 

approximate  theoretical 


analysis  yields  for  the  growth  rate. 


y~  +  2k  —  (ut,  —  ua)y  ~  kgA  =  0 

Pa  j? 

In  the  limit  where  —  »  1  and  — * »  1  this  reduces  to 
Pb  ku‘ 

y  =  4kg  -  kua 

which  is  identical  to  the  results  of  Bodner,7 


(17a) 


(17b) 


In  order  to  test  the  vorticity  generating  code,  the  case  shown  in  Fig.  D2a.l  was  run  with  pu  -  0. 
In  Fig.  D2a.2,  the  analytical  dispersion  relation  Eq.  (16a)  is  plotted  as  the  solid  line  and  our  numerical 
results  as  solid  circles.  As  can  be  seen,  the  agreement  is  excellent  over  a  wide  range  of  wavelengths. 
In  order  to  get  this  agreement,  it  was  found  that  very  fine  resolution  (Ax  =  10'6  cm)  i*  the  unstable 
region  was  needed  at  the  shortest  wavelength.  This  was  accomplished  by  utilizing  a  variably-spaced  grid 
with  the  finest  resolution  in  the  unstable  region  and  with  the  coarsest  resolution  at  the  boundaries. 

Also  shown  in  Fig.  D2a.2  are  the  results  for  pu  *•  105  cm~2  sec"1.  The  dashed  curve  is  a  plot  of 
Eq.  (1 7a)  and  the  open  circles  are  our  numerical  results.  Agreement  is  reasonable  except  at  the  shor¬ 
test  wavelengths  where  timestep  limitations  would  not  allow  good  enough  resolution. 
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Satisfied  that  our  code  was  running  correctly,  we  next  attempted  to  analyze  the  stability  of  the 
different  kinds  of  ablation  layer  profiles  that  were  discussed  previously.  This  was  done  for  a  15/a  C-H 
foil  irradiated  with  a  1  x  1013  w/cm2  beam  and  having  an  80%  absorption  rate.  We  chose  two  profiles 
of  Fig.  Dl.2  corresponding  to  an  entropy  of  4  x  1012.  The  first  profile  has  a  monotonically  rising  pres¬ 
sure  profile  and  hence  a  broad  Rayleigh-Taylor  unstable  region  (—  100/a).  The  other  pressure  profile 
peaks  just  beyond  the  density  peak  and  decreases  from  there  to  the  critical  surface.  Hence,  it  has  a 
very  narrow  unstable  region  (—  .03^).  The  vorticity  amplitudes  for  the  first  case  are  shown  in  Fig. 
D2a.3  and  for  the  second  case  in  Fig.  D2a.4.  Whereas  the  first  case  show  exponential  growth,  the 
second  case  is  stable. 


est  Dispersion  Relotions 


Fig.  D2a  2  —Comparison  ol  calculations  to  the  theoretical  Rayleigh  I  ay  lor  dispersion  relations  of  equations  (la)  (solid  line)  and 
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CASE  #2  R-T  STABLE 
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Fig.  D2a.4— Mode  amplitudes  for  the  same  external  parameters  as  Fig.  D2a.3  but  in  the  stable  region  where 

Vn  -  .04,  <V0  -  4560.  and  Sn  -  4  x  10IJ. 
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familiar  form  of  Eq.  (16a)  where 
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1  dP 
p  dx  8' 


(23) 


For  the  case  of  a  time  independent  zero  order  state  with  finite  mass  flow  Eqs.  (19a),  (I9b),  and 
(19c)  can  be  somewhat  simplified  by  use  of  the  zero  order  continuity  and  momentum  equations 


and 
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dx 
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By  subtracting  w  times  Eq.  (19b)  from  (19c)  and  substituting  for  the  gradient  of  pressure,  Eq.  (19c) 
may  be  written  as 
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and  dropping  the  primes  Eqs.  (19)  become 
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Defining 

exp  | J*  ^  tix 

and  dropping  the  primes  Eqs.  (28)  can  be  transformed  into  two  coupled  second  order  differential  eigen¬ 
value  equations. 


(29) 


d_ 

dx 


pu  exp 

d_ 
dx 


-2i L 

g 

4 

dx 


dif 

dx 


-  k2pu  exp  I  -  i  -  -  exp 


-  f  *  dx  - 

u 


n 


?_  iL  A 

g  dx 


—  -  gk 2  ^£-  exp  f^-dx  6  +  -7-  yu~^~  exp  dx  —■ 

dx  J  u  I  dx  dx  I*7  u  dx 


(30a) 


(30b) 


And  from  Eqs.  (24)  and  (25) 
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/  -  y~  +mf>u  (31) 

J  u  pug 

Equations  (30)  can  then  be  solved  numerically  for  a  given  zero-order  state. 
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PLANS  FOR  THE  EXTENSION  AND  CONTINUATION  OF  OUR  RESEARCH  IN  INERTIAL 
CONFINEMENT  FUSION 

We  plan  to  cntinue  our  research  in  three  areas.  The  first  is  a  continued  development  of  the 
Lagrangian  triangular  grid  techniques  including  compressible  flow  modifications  and  diffusive  transport. 
These  generalizations  are  logical  extensions  of  our  previous  development  work  in  the  geometry  control 
and  re-configuration  techniques  of  the  triangular  grid  algorithms.  Along  with  this  we  plan  to  investigate 
non-local  techniques  for  solving  certain  general  types  of  elliptic  problems  because  the  nature  of  our  tri¬ 
angular  grid  techniques  tends  to  defeat  most  of  the  known  fast  solution  algorithms  for  this  class  of 
problems.  The  second  area  is  to  complete  our  analysis  of  the  stability  of  shock  collapse  and  compare 
our  analytic  solutions  with  our  computational  model  based  on  the  Chester,  Chisnell,  and  Whitham 
approximation.  We  also  plan  to  extend  our  closed  form  solution  for  the  linear  stability  of  the 
"GuderleyMike  imploding  shock  into  the  nonlinear  regime  through  use  of  the  numerical  model. 
Finally,  the  third  area  is  to  continue  the  development  and  application  of  our  vorticity  generation  model 
(VGM)  to  the  problem  of  Rayleigh-Taylor  stability  in  the  laser-driven  ablation  layer.  We  plan  to  relax 
the  incompressibility  constraint  on  the  perturbation  and  examine  the  effects  of  finite  sound  speed.  We 
also  plan  to  use  the  VGM  to  help  resolve  many  of  the  ambiguities  in  previous  Rayleigh-Taylor  stability 
investigations  due  to  the  inability  to  properly  include  convection  in  the  analytic  analyses  and  the  inabil¬ 
ity  to  allow  high  enough  resolution  in  the  numerical  studies. 
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